From b068561df69161302c60e39ee7e396632fe37b8a Mon Sep 17 00:00:00 2001 From: kir Date: Tue, 28 Oct 2003 19:53:36 +0000 Subject: [PATCH] First squeeze of cleaning --- RICH/AliRICHRecon.cxx | 708 +++++------------------------------------- RICH/AliRICHRecon.h | 586 +++++++++++++++------------------- 2 files changed, 333 insertions(+), 961 deletions(-) diff --git a/RICH/AliRICHRecon.cxx b/RICH/AliRICHRecon.cxx index 6c34b54604b..8936a39973b 100644 --- a/RICH/AliRICHRecon.cxx +++ b/RICH/AliRICHRecon.cxx @@ -17,215 +17,59 @@ #include "AliRICHRecon.h" #include "AliRICHParam.h" #include +#include #include #include -#include #include -#include -#include #include #include #include #include #include #include -#include -#include -#include -#include -#include #include #include -#include #include -#include +#include #define NPointsOfRing 201 -// Geometry of the RICH at Star... static const Int_t nPadX = AliRICHParam::NpadsY(); static const Int_t nPadY = AliRICHParam::NpadsX(); static const Float_t PadSizeX = AliRICHParam::PadSizeY(); static const Float_t PadSizeY = AliRICHParam::PadSizeX(); -static const Float_t spacer = AliRICHParam::DeadZone(); -static const Float_t degree = 180/3.1415926535; -static const Float_t pi = TMath::Pi(); static const Float_t RadiatorWidth = AliRICHParam::FreonThickness(); static const Float_t QuartzWidth = AliRICHParam::QuartzThickness(); static const Float_t GapWidth = AliRICHParam::RadiatorToPads(); static const Float_t fDTheta = 0.001; // Step for sliding window -//static const Float_t fWindowWidth = 0.040; // Hough width of sliding window static const Float_t fWindowWidth = 0.060; // Hough width of sliding window -static const Int_t fThetaBin = 750; // Hough bins -static const Float_t fThetaMin = 0.0; // Theta band min -static const Float_t fThetaMax = 0.75; // Theta band max - -static const Float_t Xmin = -AliRICHParam::PcSizeY()/2.; -static const Float_t Xmax = AliRICHParam::PcSizeY()/2.; -static const Float_t Ymin = -AliRICHParam::PcSizeX()/2.; -static const Float_t Ymax = AliRICHParam::PcSizeX()/2.; - - -// Global variables... - -Bool_t fDebug = kFALSE; -Bool_t kDISPLAY = kFALSE; -Bool_t kWEIGHT = kFALSE; -Bool_t kBACKGROUND = kFALSE; -Bool_t kMINIMIZER = kFALSE; -// - -Int_t TotEvents = 0; -static Float_t xGraph[3000],yGraph[3000]; - -static Int_t NRings = 0; -static Int_t NevTOT = 0; TMinuit *gMyMinuit ; void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); - -// Float_t fEmissionPoint; -// Float_t fTrackTheta; -// Float_t fTrackPhi; -// Float_t fXtoentr; -// Float_t fYtoentr; - -// - -static TFile *outputfile; - -static TH1F *h1_photons,*h1_photacc,*h1_hough; -static TH2F *h2_tvsppos, *h2_tvspneg,*h2_func; - -static TH2F *h2_disp; - -static TH2F *h2_test2, *h2_testmap; -//static TH2F *h2_test1, *h2_test4; -static TH2F *h2_dist_p; - -static TH1F *h1_photons1, *h1_photons2; -static TH1F *h1_houghpos, *h1_houghneg; -static TH1F *h1_mass; -static TH2F *h2_mvsp; - -static TH1F *h1_hcs, *h1_hcsw; - -static TH1F *h1_nprotons; - -static TProfile *hp_1pos, *hp_1neg; -static TProfile *hp_1posnorm, *hp_1negnorm; -static TH2F *h2_1pos, *h2_1neg; -static TH2F *h2_1posnorm, *h2_1negnorm; -static TH2F *h2_mvst; - -static TH1F *h1_deltap, *h1_deltapop; -static TH1F *h1_diffTrackTheta, *h1_diffTrackPhi; -static TH1F *h1_photaccspread; - -static TH2F *h2_diffpos, *h2_diffneg; -static TH2F *h2_map, *h2_mapw; - -static TH1F *photris; - -static TH1F *gChargeMipH1, *gMultMipH1; - -static TNtuple *hn; - -static TCanvas *StarCanvas,*Display; -//static TCanvas *Displayhcs; -static TGraph *gra; -/* -static TLine *line; -static TPolyLine *poll; -static TPolyMarker *polm; -static TMarker *Point, *TrackPoints, *Photon, *PhotonAcc; -static TText *text; -*/ - -AliRICHRecon::AliRICHRecon(const char*, const char*) +//__________________________________________________________________________________________________ +AliRICHRecon::AliRICHRecon(const char*name, const char*title) + :TTask(name,title) { - + fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; + fXmin = -AliRICHParam::PcSizeY()/2.; + fXmax = AliRICHParam::PcSizeY()/2.; + fYmin = -AliRICHParam::PcSizeX()/2.; + fYmax = AliRICHParam::PcSizeX()/2.; + fRich = (AliRICH*)gAlice->GetDetector("RICH"); - InitRecon(); -} - -void AliRICHRecon::InitRecon() -{ - - outputfile = new TFile("Anal.root","RECREATE","My Analysis histos"); - if(kDISPLAY) Display = new TCanvas("Display","RICH Display",0,0,1200,750); - - h1_photons = new TH1F("h1_photons","photons",750,0.,0.75); - h1_photacc = new TH1F("h1_photacc","photons",750,0.,0.75); - h1_hough = new TH1F("h1_hough","hough",750,0.,0.75); - h1_houghpos= new TH1F("h1_houghpos","hough",750,0.,0.75); - h1_houghneg= new TH1F("h1_houghneg","hough",750,0.,0.75); - - h2_tvsppos = new TH2F("h2_tvsppos","thetac vs p",100,0.,5.,750,0.,0.75); - h2_tvspneg = new TH2F("h2_tvspneg","thetac vs p",100,0.,5.,750,0.,0.75); - h2_func = new TH2F("h2_func"," func ",800,0.,0.8,100,-100.,100.); - h2_mvsp = new TH2F("h2_mvsp","mass vs p",100,0.,5.,200,0.,2.); - h2_mvst = new TH2F("h2_mvst","mass vs t",750,0.,0.75,200,0.,2.); - h2_map = new TH2F("h2_map","h2_map",160,0.,160.,96,0.,96.); - h2_mapw = new TH2F("h2_mapw","h2_mapw",160,0.,160.,96,0.,96.); - - h2_dist_p = new TH2F("h2_dist_p","h2_dist_p",100,0.,5.,100,0.,5.); - // - - h2_disp = new TH2F("h2_disp","STAR-RICH Event Display",165,Xmin,Xmax,100,Ymin,Ymax); - - // h2_test1 = new TH2F("h2_test1","test1 map",165,-64.,64.,100,-42.,42.); - h2_test2 = new TH2F("h2_test2","test2 map",165,-64.,64.,100,-42.,42.); - // h2_test4 = new TH2F("h2_test4","test4 map",165,-64.,64.,100,-42.,42.); - h2_testmap= new TH2F("h2_testmap","test map",165,-64.,64.,100,-42.,42.); - - // - h1_photons1 = new TH1F("h1_photons1","photons",750,0.,0.75); - h1_photons2 = new TH1F("h1_photons2","photons",750,0.,0.75); - // - h1_hcs = new TH1F("h1_hcs","hcs",750,0.,750.); - h1_hcsw = new TH1F("h1_hcsw","hcsw",750,0.,750.); - // - h1_nprotons = new TH1F("h1_nprotons","n prot",30,0.,30.); - // - hp_1pos = new TProfile("hp_1pos","Nphot vs thetac pos",250,0.,0.75); - hp_1neg = new TProfile("hp_1neg","Nphot vs thetac neg",250,0.,0.75); - hp_1posnorm = new TProfile("hp_1posnorm","Nphot vs thetac pos norm",250,0.,0.75); - hp_1negnorm = new TProfile("hp_1negnorm","Nphot vs thetac neg norm",250,0.,0.75); - // - h2_1pos = new TH2F("h2_1pos","Nphot vs p pos",100,0.,5.,30,0.,30.); - h2_1neg = new TH2F("h2_1neg","Nphot vs p neg",100,0.,5.,30,0.,30.); - h2_1posnorm = new TH2F("h2_1posnorm","Nphot vs p pos norm",100,0.,5.,30,0.,30.); - h2_1negnorm = new TH2F("h2_1negnorm","Nphot vs p neg norm",100,0.,5.,30,0.,30.); - - h1_deltap = new TH1F("h1_deltap","delta_p",200,-0.5,0.5); - h1_deltapop = new TH1F("h1_deltapop","deltapop",200,-1.,1.); - h1_diffTrackTheta = new TH1F("h1_diffTrackTheta","delta theta",200,-0.25,0.25); - h1_diffTrackPhi = new TH1F("h1_diffTrackPhi","delta phi",200,-0.25,0.25); - - h1_photaccspread = new TH1F("h1_photaccspread","photons spread",200,-0.1,0.1); - - // - - h1_mass = new TH1F("h1_mass","mass",200,0.,2.); - photris = new TH1F("photris","photris",1000,0.,1.); - h2_diffneg = new TH2F("h2_diffneg","diff neg",100,-2.5,2.5,100,-2.5,2.5); - h2_diffpos = new TH2F("h2_diffpos","diff pos",100,-2.5,2.5,100,-2.5,2.5); - gChargeMipH1 = new TH1F("gChargeMipH1"," Charge Mip ",2000,0.,2000.); - gMultMipH1 = new TH1F("gMultMipH1"," Cluster Size Mip ",50,0.,50.); - - - hn = new TNtuple("hn","ntuple", + fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos"); + if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750); + fNtuple=new TNtuple("hn","ntuple", "Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge:ThetaCerenkov:NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ"); } - +//__________________________________________________________________________________________________ void AliRICHRecon::StartProcessEvent() { @@ -236,10 +80,10 @@ void AliRICHRecon::StartProcessEvent() SetFreonScaleFactor(0.994); - if(kDISPLAY) + if(fIsDISPLAY) { DrawEvent(0); -// waiting(); +// Waiting(); } Rich()->GetLoader()->LoadHits(); @@ -287,8 +131,6 @@ void AliRICHRecon::StartProcessEvent() } cout << " iPrim " << iPrim << " pHit " << pHit << endl; -// if(iPrim==0) return; -// if(iPrim>1) Fatal("StartProcessEvent"," problems with prim to hit!!! = %3i", iPrim); if (!pHit) return; @@ -337,8 +179,6 @@ void AliRICHRecon::StartProcessEvent() Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j]; Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j]; -// cout << " cluster x " << clusX[pHit->Chamber()-1][j] << " hit track x " << primLocal.X(); -// cout << " cluster y " << clusY[pHit->Chamber()-1][j] << " hit track y " << primLocal.Y() << endl; Float_t diff = sqrt(diffx*diffx + diffy*diffy); @@ -355,12 +195,6 @@ void AliRICHRecon::StartProcessEvent() cout << " diffx " << diffx << " diffy " << diffy << endl; - if(q>0) - { - h2_diffpos->Fill(diffx,diffy); - } else { - h2_diffneg->Fill(diffx,diffy); - } SetMipIndex(MaxInd); SetTrackIndex(i); @@ -380,20 +214,15 @@ void AliRICHRecon::StartProcessEvent() Int_t qch = clusQ[pHit->Chamber()-1][MaxInd]; - gChargeMipH1->Fill(qch); - gMultMipH1->Fill((Float_t)clusMul[pHit->Chamber()-1][MaxInd]); if(MinDist < 3.0 && qch > 120 && MaxInd !=0) { - if(kBACKGROUND) + if(fIsBACKGROUND) { - Float_t Xrndm = Xmin + (Xmax-Xmin)*gRandom->Rndm(280964); - Float_t Yrndm = Ymin + (Ymax-Ymin)*gRandom->Rndm(280964); - - cout << " Xrndm " << Xrndm << " Yrndm " << Yrndm << endl; - + Float_t Xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964); + Float_t Yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964); SetShiftX(Xrndm); SetShiftY(Yrndm); @@ -411,7 +240,7 @@ void AliRICHRecon::StartProcessEvent() Float_t DiffTrackTheta = 999.; Float_t DiffTrackPhi = 999.; - while( kMINIMIZER && GetHoughPhotons() > 2 + while(fIsMINIMIZER && GetHoughPhotons() > 2 && DiffNPhotons !=0 && DiffTrackTheta > 0.0001 && Nsteps < 10) @@ -426,7 +255,7 @@ void AliRICHRecon::StartProcessEvent() PatRec(); - DiffNPhotons = abs(HoughPhotonsBefore - GetHoughPhotons()); + DiffNPhotons = TMath::Abs(HoughPhotonsBefore - GetHoughPhotons()); Float_t TrackThetaAfter = GetTrackTheta(); Float_t TrackPhiAfter = GetTrackPhi(); @@ -453,38 +282,29 @@ void AliRICHRecon::StartProcessEvent() FillHistograms(); - if(kDISPLAY) DrawEvent(1); + if(fIsDISPLAY) DrawEvent(1); - waiting(); + Waiting(); } } - // - if(kDISPLAY) Display->Print("display.ps"); -} - - + if(fIsDISPLAY) fDisplay->Print("display.ps"); +}//StartProcessEvent() +//__________________________________________________________________________________________________ void AliRICHRecon::EndProcessEvent() -{ -// function called at the end of the event loop - - printf("Processed events: %d Total events: %d \n",TotEvents,NevTOT); +{// function called at the end of the event loop - outputfile->Write(); - outputfile->Close(); + fOutFile->Write(); + fOutFile->Close(); } - +//__________________________________________________________________________________________________ void AliRICHRecon::PatRec() { - cout << " in PatRec:: " << endl; Float_t TrackTheta = GetTrackTheta(); Float_t TrackPhi = GetTrackPhi(); Float_t pmod = GetTrackMomentum(); - // Int_t q = GetTrackCharge(); - - // Int_t TrackIndex = GetTrackIndex(); Int_t MipIndex = GetMipIndex(); Bool_t kPatRec = kFALSE; @@ -517,7 +337,6 @@ void AliRICHRecon::PatRec() if (j == MipIndex) continue; - // h2_test1->Fill(CandidatePhotonX[j],CandidatePhotonY[j]); if(CandidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */ @@ -554,11 +373,6 @@ void AliRICHRecon::PatRec() CandidatePhotons++; - // fill histograms - - // h2_test4->Fill(CandidatePhotonX[j],CandidatePhotonY[j]); - - // if(kDISPLAY) h1_photons->Fill(ThetaPhotonCerenkov); } @@ -573,7 +387,7 @@ void AliRICHRecon::PatRec() HoughResponse(); - NRings++; + fNrings++; FlagPhotons(); Int_t NPhotonHough = GetHoughPhotons(); @@ -585,7 +399,7 @@ void AliRICHRecon::PatRec() return; } - if(kWEIGHT) FindWeightThetaCerenkov(); + if(fIsWEIGHT) FindWeightThetaCerenkov(); Float_t ThetaCerenkov = GetThetaCerenkov(); @@ -614,9 +428,8 @@ void AliRICHRecon::PatRec() if(fDebug) { cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl; - cout << " Rings found " << NRings << " with thetac " << ThetaCerenkov << endl; + cout << " Rings found " << fNrings << " with thetac " << ThetaCerenkov << endl; - h1_hough->Fill(ThetaCerenkov,1.); cout << " Nphotons " << GetPhotonsNumber() << " Hough " << NPhotonHough @@ -643,7 +456,6 @@ void AliRICHRecon::PatRec() if(GetPhotonFlag() == 2) { - if(pmod>2.5&&ThetaCerenkov>0.65) photris->Fill(eta); xmean += eta; x2mean += eta*eta; @@ -855,8 +667,8 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t ThetaCerenkov) if(fDebug) cout << " Theta sol 1 " << ThetaSol1 << " Theta sol 2 " << ThetaSol2 << endl; - if(ThetaSol1>0 && ThetaSol1 < pi) ThetaAtQuartz = (Float_t)ThetaSol1; - if(ThetaSol2>0 && ThetaSol2 < pi) ThetaAtQuartz = (Float_t)ThetaSol2; + if(ThetaSol1>0 && ThetaSol1 < TMath::Pi()) ThetaAtQuartz = (Float_t)ThetaSol1; + if(ThetaSol2>0 && ThetaSol2 < TMath::Pi()) ThetaAtQuartz = (Float_t)ThetaSol2; SetThetaAtQuartz(ThetaAtQuartz); } @@ -1028,10 +840,6 @@ void AliRICHRecon::FindAreaAndPortionOfRing() Int_t Zone = CheckDetectorAcceptance(); -// cout << " XPointing " << XPointRing << " YPointing " << YPointRing << " Zone " << Zone << endl; -// cout << " ShiftX " << ShiftX << " ShiftY " << ShiftY << endl; -// cout << " GetXPointOnCathode() " << GetXPointOnCathode() << endl; -// cout << " GetYPointOnCathode() " << GetYPointOnCathode() << endl; if (Zone != 0) { @@ -1066,7 +874,6 @@ void AliRICHRecon::FindAreaAndPortionOfRing() Float_t PortionOfRing = ((Float_t)NPsiAccepted)/((Float_t)(NPsiTotal)); - // cout << " Area " << Area << " Portion of ring " << PortionOfRing << endl; SetAreaOfRing(Area); SetPortionOfRing(PortionOfRing); @@ -1114,36 +921,36 @@ void AliRICHRecon::FindIntersectionWithDetector() y1 = YPoint; } // - XIntersect = Xmax; + XIntersect = fXmax; YIntersect = m*(XIntersect - x0) + y0; - if (YIntersect >= Ymin && YIntersect <= Ymax && XIntersect >= x1 && XIntersect <= x2) + if (YIntersect >= fYmin && YIntersect <= fYmax && XIntersect >= x1 && XIntersect <= x2) { SetIntersectionX(XIntersect); SetIntersectionY(YIntersect); return; } // - XIntersect = Xmin; + XIntersect = fXmin; YIntersect = m*(XIntersect - x0) + y0; - if (YIntersect >= Ymin && YIntersect <= Ymax && XIntersect >= x1 && XIntersect <= x2) + if (YIntersect >= fYmin && YIntersect <= fYmax && XIntersect >= x1 && XIntersect <= x2) { SetIntersectionX(XIntersect); SetIntersectionY(YIntersect); return; } // - YIntersect = Ymax; + YIntersect = fYmax; XIntersect = (YIntersect - y0)/m + x0; - if (XIntersect >= Xmin && XIntersect <= Xmax && YIntersect >= y1 && YIntersect <= y2) + if (XIntersect >= fXmin && XIntersect <= fXmax && YIntersect >= y1 && YIntersect <= y2) { SetIntersectionX(XIntersect); SetIntersectionY(YIntersect); return; } // - YIntersect = Ymin; + YIntersect = fYmin; XIntersect = (YIntersect - y0)/m + x0; - if (XIntersect >= Xmin && XIntersect <= Xmax && YIntersect >= y1 && YIntersect <= y2) + if (XIntersect >= fXmin && XIntersect <= fXmax && YIntersect >= y1 && YIntersect <= y2) { SetIntersectionX(XIntersect); SetIntersectionY(YIntersect); @@ -1151,13 +958,9 @@ void AliRICHRecon::FindIntersectionWithDetector() } cout << " sono fuori!!!!!!" << endl; -// cout << " x1 " << x1 << " x2 " << x2 << endl; -// cout << " y1 " << y1 << " y2 " << y2 << endl; -// cout << " Xmin " << Xmin << " Xmax " << Xmax << endl; -// cout << " Ymin " << Ymin << " Ymax " << Ymax << endl; } - +//__________________________________________________________________________________________________ Int_t AliRICHRecon::CheckDetectorAcceptance() { @@ -1167,114 +970,33 @@ Int_t AliRICHRecon::CheckDetectorAcceptance() Float_t Xcoord = GetDetectorWhereX(); Float_t Ycoord = GetDetectorWhereY(); -// cout << " Xcoord " << Xcoord << " Ycoord " << Ycoord << endl; - if(Xcoord > Xmax) + if(Xcoord > fXmax) { - if(Ycoord > Ymax) return 2; - if(Ycoord > Ymin && Ycoord < Ymax) return 3; - if(Ycoord < Ymin) return 4; + if(Ycoord > fYmax) return 2; + if(Ycoord > fYmin && Ycoord < fYmax) return 3; + if(Ycoord < fYmin) return 4; } - if(Xcoord < Xmin) + if(Xcoord < fXmin) { - if(Ycoord > Ymax) return 8; - if(Ycoord > Ymin && Ycoord < Ymax) return 7; - if(Ycoord < Ymin) return 6; + if(Ycoord > fYmax) return 8; + if(Ycoord > fYmin && Ycoord < fYmax) return 7; + if(Ycoord < fYmin) return 6; } - if(Xcoord > Xmin && Xcoord < Xmax) + if(Xcoord > fXmin && Xcoord < fXmax) { - if(Ycoord > Ymax) return 1; - if(Ycoord > Ymin && Ycoord < Ymax) return 0; - if(Ycoord < Ymin) return 5; + if(Ycoord > fYmax) return 1; + if(Ycoord > fYmin && Ycoord < fYmax) return 0; + if(Ycoord < fYmin) return 5; } return 999; } - -void AliRICHRecon::DrawRing() -{ - - // Float_t xGraph[1000],yGraph[1000]; - - Float_t type; - // Float_t MassOfParticle; - Float_t beta; - Float_t nfreon; - - Float_t ThetaCerenkov; - - // Float_t Xtoentr = GetEntranceX(); - // Float_t Ytoentr = GetEntranceY(); - - // Float_t pmod = GetTrackMomentum(); - // Float_t TrackTheta = GetTrackTheta(); - // Float_t TrackPhi = GetTrackPhi(); - - SetPhotonEnergy(6.85); - SetFreonRefractiveIndex(); - - SetEmissionPoint(RadiatorWidth/2.); - - type = 1; - - if(type == 1) - { - SetMassHypotesis(0.139567); - SetBetaOfParticle(); - - beta = GetBetaOfParticle(); - - } - else if(type == 2) - { - ThetaCerenkov = GetThetaCerenkov(); - FindBetaFromTheta(ThetaCerenkov); - } - - nfreon = GetFreonRefractiveIndex(); - - Float_t thetacer = Cerenkovangle(nfreon,beta); - - if(fDebug) cout << " TetaCer in DrawRing " << thetacer << endl; - - Int_t nPoints = 100; - - Int_t nPointsToDraw = 0; - for(Int_t i=0;iDraw("same"); - gra = new TGraph(nPointsToDraw,xGraph,yGraph); - gra->Draw("AC"); - StarCanvas->Update(); - -} - +//__________________________________________________________________________________________________ Float_t AliRICHRecon::PhotonPositionOnCathode() { // Float_t MassOfParticle; Float_t beta; Float_t nfreon; - // Float_t pmod = GetTrackMomentum(); - // Float_t TrackTheta = GetTrackTheta(); - // Float_t TrackPhi = GetTrackPhi(); - - // Float_t phpad = GetPhiPoint(); SetPhotonEnergy(6.85); SetEmissionPoint(RadiatorWidth/2.); @@ -1286,17 +1008,10 @@ Float_t AliRICHRecon::PhotonPositionOnCathode() beta = GetBetaOfParticle(); nfreon = GetFreonRefractiveIndex(); - // Float_t thetacer = Cerenkovangle(nfreon,beta); - - // cout << " FromEmissionToCathode: thetacer " << thetacer << " phpad " << phpad << endl; Float_t Radius = FromEmissionToCathode(); if (Radius == 999.) return 999.; - // Float_t Xphoton = GetXPointOnCathode(); - // Float_t Yphoton = GetYPointOnCathode(); - // cout << " PhotonPositionOnCathode: Xphoton " << Xphoton << " Yphoton " << Yphoton << - // " Radius for photon " << Radius << endl; return 0; } @@ -1386,8 +1101,6 @@ Float_t AliRICHRecon::FromEmissionToCathode() Float_t yq = QuartzWidth*sin(Phi)*tan(thetaquar); Float_t yg = GapWidth*sin(Phi)*tan(thetagap); -// Float_t xtot = x1 + xw + xq + xg; -// Float_t ytot = y1 + yw + yq + yg; Float_t xtot = xEmiss + xw + xq + xg; Float_t ytot = yEmiss + yw + yq + yg; @@ -1395,7 +1108,6 @@ Float_t AliRICHRecon::FromEmissionToCathode() SetXPointOnCathode(xtot); SetYPointOnCathode(ytot); -// cout << " xtot " << xtot << " ytot " << ytot << endl; Float_t DistanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2) +TMath::Power(fPhotonLimitY,2)); @@ -1517,7 +1229,7 @@ void AliRICHRecon::HoughResponse() // calculate weights - if(kWEIGHT) + if(fIsWEIGHT) { lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta; SetThetaOfRing(lowerlimit); @@ -1555,8 +1267,6 @@ void AliRICHRecon::HoughResponse() // cout << "weight..." << weight << endl; - h1_photons1->Fill(angle); - h1_photons2->Fill(angle,weight); if (bin1<0) bin1=0; if (bin2>nBin) bin2=nBin; @@ -1569,14 +1279,6 @@ void AliRICHRecon::HoughResponse() } } -// if(kDISPLAY) -// { -// for(Int_t j=0;j<750;j++) -// { -// h1_hcs->Fill(((Float_t)j),hcs[j]); -// h1_hcsw->Fill(((Float_t)j),hcsw[j]); -// } -// } if(WeightFlag == 0) { @@ -1728,198 +1430,6 @@ void AliRICHRecon::DrawEvent(Int_t flag) { flag=1; // dummy to be removed... -/* - Float_t xGraph[3000],yGraph[3000]; - - Float_t ThetaCerenkov; - - // Display event... - - gStyle->SetPalette(1,0); - - if(flag == 0) - { - - // Display = new TCanvas("Display","Star Display",0,0,1200,750); - - Display->ToggleEventStatus(); - Display->Modified() - - text = new TText(0,0,""); - text->SetTextFont(61); - text->SetTextSize(0.03); - text->SetTextAlign(22); - - Display->Resize(); - - h2_disp->Reset(); - - for(Int_t j=1;j<=nPixels;j++) - { - Float_t xpad = fPixels_localX[j-1]; - Float_t ypad = fPixels_localY[j-1]; - h2_disp->Fill(xpad,ypad,fPixels_charge[j-1]); - } - - h2_disp->SetMaximum(200); - // h2_disp->SetMaximum(1); - h2_disp->SetStats(0); - h2_disp->Draw("colz"); - - for(Int_t i=0; iSetMarkerSize(1.5); - - Float_t pmod = sqrt(fRichPrimaries_localPadPx[i] * fRichPrimaries_localPadPx[i] + - fRichPrimaries_localPadPy[i] * fRichPrimaries_localPadPy[i] + - fRichPrimaries_localPadPz[i] * fRichPrimaries_localPadPz[i]); - - if(pmod < 1) TrackPoints->SetMarkerColor(kBlue); - if(pmod > 1 && pmod < 2) TrackPoints->SetMarkerColor(kGreen); - if(pmod > 2) TrackPoints->SetMarkerColor(kRed); - - TrackPoints->Draw(); - - line = new TLine(-0.13,-42.,-0.13,42.); - line->Draw(); - line = new TLine(0.13,-42.,0.13,42.); - line->Draw(); - line = new TLine(-64.,-0.13,64.,-0.13); - line->Draw(); - line = new TLine(-64.,0.13,64.,0.13); - line->Draw(); - - } - - return; - - } - - // - - // Draw rings... - - // - - // Float_t Xtoentr = GetEntranceX(); - // Float_t Ytoentr = GetEntranceY(); - - // Float_t pmod = GetTrackMomentum(); - // Float_t TrackTheta = GetTrackTheta(); - // Float_t TrackPhi = GetTrackPhi(); - - SetPhotonEnergy(6.85); - SetFreonRefractiveIndex(); - - SetEmissionPoint(RadiatorWidth/2.); - - ThetaCerenkov = GetThetaCerenkov(); - - if (ThetaCerenkov == 999.) return; - - Int_t nPointsToDraw = 0; - - for(Int_t i=0;i<99;i++) - { - Float_t phpad = 2*TMath::Pi()*i/99; - SetThetaPhotonInTRS(ThetaCerenkov); - SetPhiPhotonInTRS(phpad); - FindPhotonAnglesInDRS(); - Float_t Radius = FromEmissionToCathode(); - - if (Radius == 999.) continue; - - Float_t ShiftX = GetShiftX(); - Float_t ShiftY = GetShiftY(); - - Float_t XPointRing = GetXPointOnCathode() + ShiftX; - Float_t YPointRing = GetYPointOnCathode() + ShiftY; - - SetDetectorWhereX(XPointRing); - SetDetectorWhereY(YPointRing); - - Int_t Zone = CheckDetectorAcceptance(); - - if (Zone != 0) - { - FindIntersectionWithDetector(); - xGraph[nPointsToDraw] = GetIntersectionX(); - yGraph[nPointsToDraw] = GetIntersectionY(); - nPointsToDraw++; - } - else - { - xGraph[nPointsToDraw] = GetXPointOnCathode() + GetShiftX(); - yGraph[nPointsToDraw] = GetYPointOnCathode() + GetShiftY(); - nPointsToDraw++; - } - } - - xGraph[nPointsToDraw] = xGraph[0]; - yGraph[nPointsToDraw] = yGraph[0]; - - poll = new TPolyLine(nPointsToDraw+1,xGraph,yGraph); - poll->SetLineColor(2); - poll->SetLineWidth(3); - - Display->Update(); - - // waiting(); - poll->Draw(); - - for(Int_t j=0;j= 1) -// { - -// Photon = new TMarker(xhit,yhit,4); -// Photon->SetMarkerSize(1.5); -// Photon->Draw("same"); - -// } - - - if(FlagPhoton == 2) - { - - PhotonAcc = new TMarker(xhit,yhit,30); - PhotonAcc->SetMarkerSize(1.5); - PhotonAcc->SetMarkerColor(50); - PhotonAcc->Draw("same"); - - } - } - - Display->Update(); - -// waiting(); -// h1_photons->Draw(); -// Display->Update(); - -// waiting(); -// h1_photacc->Draw(); -// Display->Update(); - -// waiting(); - -// Display->Update(); - -// h1_photons->Reset(); -// h1_photacc->Reset(); - -*/ } Float_t AliRICHRecon::FindMassOfParticle() @@ -1971,11 +1481,11 @@ void AliRICHRecon::FillHistograms() if(fDebug) { cout << " p " << pmod << " ThetaC " << ThetaCerenkov - << " rings " << NRings << endl; + << " rings " << fNrings << endl; } Int_t NPhotonHough = GetHoughPhotons(); - Float_t NPhotonHoughNorm = GetHoughPhotonsNorm(); +// Float_t NPhotonHoughNorm = GetHoughPhotonsNorm(); Float_t InRing = GetPortionOfRing(); Float_t MassOfParticle = FindMassOfParticle(); @@ -1983,12 +1493,9 @@ void AliRICHRecon::FillHistograms() Float_t HoughArea = GetHoughArea(); Float_t Multiplicity = GetEventMultiplicity(); -// cout << " area " << HoughArea << " mult " << Multiplicity << endl; Float_t var[20]; -// var[0] = (Float_t)runID; -// var[1] = (Float_t)evID; var[0] = 0; var[1] = 0; var[2] = VertZ; @@ -2010,77 +1517,25 @@ void AliRICHRecon::FillHistograms() var[18] = TPCLastZ; var[19] = MinDist; - hn->Fill(var); + fNtuple->Fill(var); - h1_mass->Fill(MassOfParticle); - h2_mvsp->Fill(pmod,MassOfParticle); - h2_mvst->Fill(ThetaCerenkov,MassOfParticle); FittedTrackTheta = GetFittedTrackTheta(); FittedTrackPhi = GetFittedTrackPhi(); - Float_t DiffTheta = FittedTrackTheta - TrackTheta; - Float_t DiffPhi = FittedTrackPhi - TrackPhi; - h1_diffTrackTheta -> Fill(DiffTheta); - h1_diffTrackPhi -> Fill(DiffPhi); - if(ThetaCerenkov > 0.505 && ThetaCerenkov < 0.605) - { + if(ThetaCerenkov > 0.505 && ThetaCerenkov < 0.605) { SetPhotonEnergy(6.85); SetFreonRefractiveIndex(); - - Float_t pmom = GetTrackMomentum(); - Float_t beta = 1./(cos(ThetaCerenkov)*GetFreonRefractiveIndex()); - Float_t gamma = 1./sqrt(1.-beta*beta); - - Float_t pmomnew = 0.93828*beta*gamma; - Float_t deltap = pmomnew - pmom; - h1_deltap->Fill(deltap); - Float_t deltapop = deltap/pmom; - h1_deltapop->Fill(deltapop); - - h1_nprotons->Fill((Float_t)NPhotonHoughNorm); - } - - if(q > 0) - { - h2_tvsppos->Fill(pmod,ThetaCerenkov); - hp_1pos->Fill(ThetaCerenkov,(Float_t)NPhotonHough); - hp_1posnorm->Fill(ThetaCerenkov,(Float_t)NPhotonHoughNorm); - h2_1pos->Fill(pmod,(Float_t)NPhotonHough); - h2_1posnorm->Fill(pmod,(Float_t)NPhotonHoughNorm); - h1_houghpos->Fill(ThetaCerenkov); - } -else - { - h2_tvspneg->Fill(pmod,ThetaCerenkov); - hp_1neg->Fill(ThetaCerenkov,(Float_t)NPhotonHough); - hp_1negnorm->Fill(ThetaCerenkov,(Float_t)NPhotonHoughNorm); - h2_1neg->Fill(pmod,(Float_t)NPhotonHough); - h2_1negnorm->Fill(pmod,(Float_t)NPhotonHoughNorm); - h1_houghneg->Fill(ThetaCerenkov); } Int_t NPhotons = GetPhotonsNumber(); for (Int_t j=0; j < NPhotons;j++) - - { - SetPhotonIndex(j); - - Float_t eta = GetPhotonEta(); - - if(GetPhotonFlag() == 2) - { - h1_photacc->Fill(eta); - Float_t photaccspread = eta - ThetaCerenkov; - h1_photaccspread->Fill(photaccspread); - } - - } -} - + SetPhotonIndex(j); +}//FillHistograms() +//__________________________________________________________________________________________________ void AliRICHRecon::Minimization() { @@ -2223,7 +1678,7 @@ void AliRICHRecon::EstimationOfTheta() SetEstimationOfThetaRMS(RMS); } -void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t iflag) +void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t) { AliRICHRecon *gMyRecon = (AliRICHRecon*)gMyMinuit->GetObjectFit(); @@ -2241,21 +1696,21 @@ void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, In f = (Double_t)(1000*RMS/(Float_t)HoughPhotons); - if(fDebug) cout << " f " << f - << " theta " << par[0] << " phi " << par[1] - << " HoughPhotons " << HoughPhotons << endl; - - if(fDebug&&iflag == 3) - { - cout << " --- end convergence...summary --- " << endl; - cout << " theta " << par[0] << endl; - cout << " phi " << par[1] << endl; - } +// if(fDebug) cout << " f " << f +// << " theta " << par[0] << " phi " << par[1] +// << " HoughPhotons " << HoughPhotons << endl; +// +// if(fDebug&&iflag == 3) +// { +// cout << " --- end convergence...summary --- " << endl; +// cout << " theta " << par[0] << endl; +// cout << " phi " << par[1] << endl; +// } } -void AliRICHRecon::waiting() +void AliRICHRecon::Waiting() { - if(!kDISPLAY) return; + if(!fIsDISPLAY) return; cout << " Press any key to continue..."; // gSystem->ProcessEvents(); @@ -2266,8 +1721,3 @@ void AliRICHRecon::waiting() return; } -/* -void ~AliRICHRecon() -{ -} -*/ diff --git a/RICH/AliRICHRecon.h b/RICH/AliRICHRecon.h index f08004aa4bf..c48ba4cf23a 100644 --- a/RICH/AliRICHRecon.h +++ b/RICH/AliRICHRecon.h @@ -5,338 +5,182 @@ * See cxx source for full Copyright notice */ #include -#include "AliRICH.h" -#include class AliRICH; - -class AliRICHRecon : public TTask { - public : - - AliRICHRecon(const char*, const char*); +class TFile; +class TNtuple; +class TCanvas; + +class AliRICHRecon : public TTask +{ +public : + AliRICHRecon(const char*, const char*); ~AliRICHRecon(){EndProcessEvent();} - AliRICH* Rich() {return fRich;} - void StartProcessEvent(); - void EndProcessEvent(); - - void InitRecon(); - - void PatRec(); - - void Minimization(); - - void FillHistograms(); - - void FindThetaPhotonCerenkov(); - void FindAreaAndPortionOfRing(); - - void FindEmissionPoint(); - - Int_t PhotonInBand(); - Float_t PhotonPositionOnCathode(); - - void FindPhotonAnglesInDRS(); - void FindPhiPoint(); - void FindThetaAtQuartz(Float_t ThetaCerenkov); - - Float_t FindMassOfParticle(); - - Float_t Cerenkovangle(Float_t n, Float_t beta); - - void HoughResponse(); - void HoughFiltering(float HCS[]); - void FlagPhotons(); - void FindWeightThetaCerenkov(); - - void EstimationOfTheta(); - void FindIntersectionWithDetector(); - Int_t CheckDetectorAcceptance(); - - void DrawRing(); - void DrawEvent(Int_t flag); - - void waiting(); - - ////////////////////////////////////// - Float_t GetEventVertexZ() { return fEventVertZ;} - Float_t GetEventMultiplicity() { return fEventMultiplicity;} - Float_t GetPhotonEnergy() { return fPhotonEnergy;} - Float_t GetFreonRefractiveIndex() { return fFreonRefractiveIndex;} - Float_t GetQuartzRefractiveIndex() { return fQuartzRefractiveIndex;} - Float_t GetGasRefractiveIndex() { return fGasRefractiveIndex;} - - Float_t GetFreonScaleFactor() { return fFreonScaleFactor;} - - Float_t GetEmissionPoint() { return fLengthEmissionPoint;} - Float_t GetMassHypotesis() { return fMassHypotesis;} - Float_t GetBetaOfParticle() { return fTrackBeta;} - Float_t GetEntranceX() { return fXtoentr;} - Float_t GetEntranceY() { return fYtoentr;} - Float_t GetThetaCerenkov() { return fThetaCerenkov;} - Float_t GetThetaPhotonCerenkov() { return fThetaPhotonCerenkov;} - Float_t GetTrackMomentum() { return fTrackMomentum;} - Float_t GetTrackEta() { return fTrackEta;} - Float_t GetTrackTheta() { return fTrackTheta;} - Float_t GetTrackPhi() { return fTrackPhi;} - Float_t GetTrackPt() { return fTrackPt;} - Int_t GetTrackCharge() { return fTrackCharge;} - Float_t GetTrackTPCLastZ() { return fTrackTPCLastZ;} - Float_t GetMinDist() { return fMinDist;} - Float_t GetXPointOnCathode() { return fPhotonLimitX;} - Float_t GetYPointOnCathode() { return fPhotonLimitY;} - - - Float_t GetThetaPhotonInDRS() { return fThetaPhotonInDRS;} - Float_t GetPhiPhotonInDRS() { return fPhiPhotonInDRS;} - Float_t GetThetaPhotonInTRS() { return fThetaPhotonInTRS;} - Float_t GetPhiPhotonInTRS() { return fPhiPhotonInTRS;} - - Float_t GetThetaAtQuartz() { return fThetaAtQuartz;} - - Float_t GetPhiPoint(){ return fPhiPoint;} - Float_t GetXCoordOfEmission() {return fXEmiss;} - Float_t GetYCoordOfEmission() {return fYEmiss;} - - Float_t GetXInnerRing() {return fXInner;} - Float_t GetYInnerRing() {return fYInner;} - Float_t GetRadiusInnerRing() {return fInnerRadius;} - - Float_t GetXOuterRing() {return fXOuter;} - Float_t GetYOuterRing() {return fYOuter;} - Float_t GetRadiusOuterRing() {return fOuterRadius;} - Float_t GetShiftX() { return fShiftX;} - Float_t GetShiftY() { return fShiftY;} - Float_t GetDetectorWhereX() { return fXcoord;} - Float_t GetDetectorWhereY() { return fYcoord;} - Float_t GetIntersectionX() {return fIntersectionX;} - Float_t GetIntersectionY() {return fIntersectionY;} - - Float_t GetThetaOfRing() { return fThetaOfRing;} - Float_t GetAreaOfRing() {return fAreaOfRing;} - Float_t GetPortionOfRing() {return fPortionOfRing;} - Float_t GetHoughArea() { return fHoughArea;} - - Int_t GetPhotonsNumber() { return fPhotonsNumber;} - Int_t GetPhotonIndex() { return fPhotonIndex;} - Float_t GetPhotonEta() { return fPhotonEta[fPhotonIndex];} - Int_t GetPhotonFlag() { return fPhotonFlag[fPhotonIndex];} - Float_t GetPhotonWeight() { return fPhotonWeight[fPhotonIndex];} - - Float_t GetHoughRMS() { return fHoughRMS;} - - Int_t GetMipIndex() { return fMipIndex;} - Int_t GetTrackIndex() { return fTrackIndex;} - Float_t* GetCandidatePhotonX() { return fCandidatePhotonX;} - Float_t* GetCandidatePhotonY() { return fCandidatePhotonY;} - Int_t GetCandidatePhotonsNumber() { return fCandidatePhotonsNumber;} - - Int_t GetHoughPhotons() { return fHoughPhotons;} - Float_t GetHoughPhotonsNorm() { return fHoughPhotonsNorm;} - - Float_t GetFittedTrackTheta() { return fFittedTrackTheta;} - Float_t GetFittedTrackPhi() { return fFittedTrackPhi;} - Float_t GetFittedThetaCerenkov() { return fFittedThetaCerenkov;} - Int_t GetFittedHoughPhotons() { return fFittedHoughPhotons;} - Float_t GetEstimationOfTheta() { return fEstimationOfTheta;} - Float_t GetEstimationOfThetaRMS() { return fEstimationOfThetaRMS;} - - void SetEventVertexZ(Float_t EventVertZ) { fEventVertZ = EventVertZ;} - void SetEventMultiplicity(Float_t EventMultiplicity) { fEventMultiplicity = EventMultiplicity;} - void SetPhotonEnergy(Float_t PhotonEnergy) { fPhotonEnergy = PhotonEnergy;} - void SetFreonRefractiveIndex() {fFreonRefractiveIndex = fFreonScaleFactor*(1.177+0.0172*fPhotonEnergy);} - void SetQuartzRefractiveIndex() {fQuartzRefractiveIndex = sqrt(1+(46.411/(113.763556-TMath::Power(fPhotonEnergy,2)))+(228.71/(328.51563-TMath::Power(fPhotonEnergy,2))));} - void SetGasRefractiveIndex() { fGasRefractiveIndex = 1.;} - - void SetFreonScaleFactor(Float_t FreonScaleFactor) {fFreonScaleFactor = FreonScaleFactor;} - - void SetEmissionPoint(Float_t LengthEmissionPoint) { fLengthEmissionPoint = LengthEmissionPoint;} - void SetMassHypotesis(Float_t mass) {fMassHypotesis = mass;} - - void SetBetaOfParticle() { fTrackBeta = fTrackMomentum/sqrt(TMath::Power(fTrackMomentum,2)+TMath::Power(fMassHypotesis,2));} - - void SetEntranceX(Float_t Xtoentr) { fXtoentr = Xtoentr;} - void SetEntranceY(Float_t Ytoentr) { fYtoentr = Ytoentr;} - - void SetThetaPhotonInTRS(Float_t Theta) {fThetaPhotonInTRS = Theta;} - void SetPhiPhotonInTRS(Float_t Phi) {fPhiPhotonInTRS = Phi;} - void SetThetaPhotonInDRS(Float_t Theta) {fThetaPhotonInDRS = Theta;} - void SetPhiPhotonInDRS(Float_t Phi) {fPhiPhotonInDRS = Phi;} - - void SetThetaAtQuartz(Float_t ThetaAtQuartz) {fThetaAtQuartz = ThetaAtQuartz;} - - void SetPhiPoint(Float_t PhiPoint){ fPhiPoint = PhiPoint;} - - void SetXCoordOfEmission(Float_t XEmiss) {fXEmiss = XEmiss;} - void SetYCoordOfEmission(Float_t YEmiss) {fYEmiss = YEmiss;} - - - void SetXPointOnCathode(Float_t PhotonLimitX) { fPhotonLimitX = PhotonLimitX;} - void SetYPointOnCathode(Float_t PhotonLimitY) { fPhotonLimitY = PhotonLimitY;} - - void SetXInnerRing(Float_t XInner) {fXInner = XInner;} - void SetYInnerRing(Float_t YInner) {fYInner = YInner;} - void SetRadiusInnerRing(Float_t InnerRadius) {fInnerRadius = InnerRadius;} - - void SetXOuterRing(Float_t XOuter) {fXOuter = XOuter;} - void SetYOuterRing(Float_t YOuter) {fYOuter = YOuter;} - void SetRadiusOuterRing(Float_t OuterRadius) {fOuterRadius = OuterRadius;} - - void SetThetaCerenkov(Float_t ThetaCer) {fThetaCerenkov = ThetaCer;} - void SetThetaPhotonCerenkov(Float_t ThetaPhotCer) {fThetaPhotonCerenkov = ThetaPhotCer;} - - void SetTrackMomentum(Float_t TrackMomentum) {fTrackMomentum = TrackMomentum;} - void SetTrackEta(Float_t TrackEta) {fTrackEta = TrackEta;} - void SetTrackTheta(Float_t TrackTheta) { fTrackTheta = TrackTheta;} - void SetTrackPhi(Float_t TrackPhi) { fTrackPhi = TrackPhi;} - void SetTrackPt(Float_t TrackPt) { fTrackPt = TrackPt;} - void SetTrackCharge(Int_t TrackCharge) { fTrackCharge = TrackCharge;} - void SetTrackTPCLastZ(Float_t TrackTPCLastZ) { fTrackTPCLastZ = TrackTPCLastZ;} - void SetMinDist(Float_t MinDist) { fMinDist = MinDist;} - void SetShiftX(Float_t ShiftX) { fShiftX = ShiftX;} - void SetShiftY(Float_t ShiftY) { fShiftY = ShiftY;} - - void SetDetectorWhereX(Float_t Xcoord) { fXcoord = Xcoord;} - void SetDetectorWhereY(Float_t Ycoord) { fYcoord = Ycoord;} - - void SetIntersectionX(Float_t IntersectionX) { fIntersectionX = IntersectionX;} - void SetIntersectionY(Float_t IntersectionY) { fIntersectionY = IntersectionY;} - - void SetThetaOfRing(Float_t ThetaOfRing) { fThetaOfRing = ThetaOfRing;} - void SetAreaOfRing(Float_t AreaOfRing) { fAreaOfRing = AreaOfRing;} - void SetPortionOfRing(Float_t PortionOfRing) { fPortionOfRing = PortionOfRing;} - void SetHoughArea(Float_t HoughArea) { fHoughArea = HoughArea;} - - - void SetPhotonsNumber(Int_t PhotonsNumber) { fPhotonsNumber = PhotonsNumber;} - void SetPhotonIndex(Int_t PhotonIndex) { fPhotonIndex = PhotonIndex;} - void SetPhotonEta(Float_t PhotonEta) { fPhotonEta[fPhotonIndex] = PhotonEta;} - void SetPhotonFlag(Int_t PhotonFlag) { fPhotonFlag[fPhotonIndex] = PhotonFlag;} - void SetPhotonWeight(Float_t PhotonWeight) { fPhotonWeight[fPhotonIndex] = PhotonWeight;} - - void SetHoughRMS(Float_t HoughRMS) { fHoughRMS = HoughRMS;} - void SetMipIndex(Int_t MipIndex) { fMipIndex = MipIndex;} - void SetTrackIndex(Int_t TrackIndex) { fTrackIndex = TrackIndex;} - - void SetCandidatePhotonX(Float_t *CandidatePhotonX) { fCandidatePhotonX = CandidatePhotonX;} - void SetCandidatePhotonY(Float_t *CandidatePhotonY) { fCandidatePhotonY = CandidatePhotonY;} - void SetCandidatePhotonsNumber(Int_t CandidatePhotonsNumber) { fCandidatePhotonsNumber = CandidatePhotonsNumber;} - void SetHoughPhotons(Int_t HoughPhotons) { fHoughPhotons = HoughPhotons;} - void SetHoughPhotonsNorm(Float_t HoughPhotonsNorm) { fHoughPhotonsNorm = HoughPhotonsNorm;} - - void SetFittedTrackTheta(Float_t FittedTrackTheta) { fFittedTrackTheta = FittedTrackTheta;} - void SetFittedTrackPhi(Float_t FittedTrackPhi) { fFittedTrackPhi = FittedTrackPhi;} - void SetFittedThetaCerenkov(Float_t FittedThetaCerenkov) { fFittedThetaCerenkov = FittedThetaCerenkov;} - void SetFittedHoughPhotons(Int_t FittedHoughPhotons) { fFittedHoughPhotons = FittedHoughPhotons;} - void SetEstimationOfTheta(Float_t EstimationOfTheta) { fEstimationOfTheta = EstimationOfTheta;} - void SetEstimationOfThetaRMS(Float_t EstimationOfThetaRMS) { fEstimationOfThetaRMS = EstimationOfThetaRMS;} - - void FindBetaFromTheta(Float_t ThetaCerenkov) {fTrackBeta = 1/(fFreonRefractiveIndex*cos(ThetaCerenkov));} - - Float_t SnellAngle(Float_t n1, Float_t n2, Float_t theta1); - - Float_t FromEmissionToCathode(); - ////////////////////////////////////// - private: - - AliRICH* fRich; - - Float_t fEventVertZ; - Float_t fEventMultiplicity; - - Float_t fTrackTheta; - Float_t fTrackPhi; - Float_t fTrackMomentum; - Float_t fTrackEta; - Float_t fTrackPt; - Int_t fTrackCharge; - Float_t fTrackTPCLastZ; - Float_t fMinDist; - Float_t fTrackBeta; - - Float_t fXtoentr; - Float_t fYtoentr; - - Float_t fThetaPhotonInTRS; - Float_t fPhiPhotonInTRS; - - Float_t fThetaPhotonInDRS; - Float_t fPhiPhotonInDRS; - - Float_t fThetaAtQuartz; - Float_t fPhiPoint; - - Float_t fXEmiss; - Float_t fYEmiss; - - Float_t fXInner; - Float_t fYInner; - Float_t fXOuter; - Float_t fYOuter; - Float_t fInnerRadius; - Float_t fOuterRadius; - - Float_t fPhotonEnergy; - Float_t fFreonRefractiveIndex; - Float_t fQuartzRefractiveIndex; - Float_t fGasRefractiveIndex; - - Float_t fFreonScaleFactor; - - Float_t fLengthEmissionPoint; - - Float_t fPhotonLimitX; - Float_t fPhotonLimitY; - Float_t fDistanceFromCluster; - - Float_t fMassHypotesis; - - Float_t fCerenkovAnglePad; - - Float_t fThetaPhotonCerenkov; - - Float_t fShiftX; - Float_t fShiftY; - - Float_t fXcoord; - Float_t fYcoord; - - Float_t fIntersectionX; - Float_t fIntersectionY; - - Float_t fThetaOfRing; - Float_t fAreaOfRing; - Float_t fPortionOfRing; - Float_t fHoughArea; - - Int_t fPhotonsNumber; - Int_t fPhotonIndex; - Float_t fPhotonEta[3000]; - Int_t fPhotonFlag[3000]; - Float_t fPhotonWeight[3000]; - - Float_t fHoughRMS; - - Int_t fMipIndex; - Int_t fTrackIndex; - - Float_t* fCandidatePhotonX; - Float_t* fCandidatePhotonY; - Int_t fCandidatePhotonsNumber; - - Int_t fHoughPhotons; - Float_t fHoughPhotonsNorm; - - Float_t fFittedTrackTheta; - Float_t fFittedTrackPhi; - Float_t fFittedThetaCerenkov; - Int_t fFittedHoughPhotons; - - Float_t fEstimationOfTheta; - Float_t fEstimationOfThetaRMS; - - public: - + AliRICH* Rich() {return fRich;} + void StartProcessEvent(); + void EndProcessEvent(); + void InitRecon(); + void PatRec(); + void Minimization(); + void FillHistograms(); + void FindThetaPhotonCerenkov(); + void FindAreaAndPortionOfRing(); + void FindEmissionPoint(); + void FindPhotonAnglesInDRS(); + void FindPhiPoint(); + void FindThetaAtQuartz(Float_t ThetaCerenkov); + void HoughResponse(); + void HoughFiltering(float HCS[]); + void FlagPhotons(); + void FindWeightThetaCerenkov(); + void EstimationOfTheta(); + void FindIntersectionWithDetector(); + void DrawEvent(Int_t flag); + void Waiting(); + Int_t GetPhotonFlag() { return fPhotonFlag[fPhotonIndex];} + Int_t PhotonInBand(); + Int_t CheckDetectorAcceptance(); + Int_t GetTrackCharge() { return fTrackCharge;} + Int_t GetPhotonsNumber() { return fPhotonsNumber;} + Int_t GetPhotonIndex() { return fPhotonIndex;} + Int_t GetMipIndex() { return fMipIndex;} + Int_t GetTrackIndex() { return fTrackIndex;} + Int_t GetCandidatePhotonsNumber() { return fCandidatePhotonsNumber;} + Int_t GetHoughPhotons() { return fHoughPhotons;} + Int_t GetFittedHoughPhotons() { return fFittedHoughPhotons;} + Float_t FindMassOfParticle(); + Float_t Cerenkovangle(Float_t n, Float_t beta); + Float_t PhotonPositionOnCathode(); + Float_t GetEventVertexZ() { return fEventVertZ;} + Float_t GetEventMultiplicity() { return fEventMultiplicity;} + Float_t GetPhotonEnergy() { return fPhotonEnergy;} + Float_t GetFreonRefractiveIndex() { return fFreonRefractiveIndex;} + Float_t GetQuartzRefractiveIndex() { return fQuartzRefractiveIndex;} + Float_t GetGasRefractiveIndex() { return fGasRefractiveIndex;} + Float_t GetFreonScaleFactor() { return fFreonScaleFactor;} + Float_t GetEmissionPoint() { return fLengthEmissionPoint;} + Float_t GetMassHypotesis() { return fMassHypotesis;} + Float_t GetBetaOfParticle() { return fTrackBeta;} + Float_t GetEntranceX() { return fXtoentr;} + Float_t GetEntranceY() { return fYtoentr;} + Float_t GetThetaCerenkov() { return fThetaCerenkov;} + Float_t GetThetaPhotonCerenkov() { return fThetaPhotonCerenkov;} + Float_t GetTrackMomentum() { return fTrackMomentum;} + Float_t GetTrackEta() { return fTrackEta;} + Float_t GetTrackTheta() { return fTrackTheta;} + Float_t GetTrackPhi() { return fTrackPhi;} + Float_t GetTrackPt() { return fTrackPt;} + Float_t GetTrackTPCLastZ() { return fTrackTPCLastZ;} + Float_t GetMinDist() { return fMinDist;} + Float_t GetXPointOnCathode() { return fPhotonLimitX;} + Float_t GetYPointOnCathode() { return fPhotonLimitY;} + Float_t GetThetaPhotonInDRS() { return fThetaPhotonInDRS;} + Float_t GetPhiPhotonInDRS() { return fPhiPhotonInDRS;} + Float_t GetThetaPhotonInTRS() { return fThetaPhotonInTRS;} + Float_t GetPhiPhotonInTRS() { return fPhiPhotonInTRS;} + Float_t GetThetaAtQuartz() { return fThetaAtQuartz;} + Float_t GetPhiPoint(){ return fPhiPoint;} + Float_t GetXCoordOfEmission() {return fXEmiss;} + Float_t GetYCoordOfEmission() {return fYEmiss;} + Float_t GetXInnerRing() {return fXInner;} + Float_t GetYInnerRing() {return fYInner;} + Float_t GetRadiusInnerRing() {return fInnerRadius;} + Float_t GetXOuterRing() {return fXOuter;} + Float_t GetYOuterRing() {return fYOuter;} + Float_t GetRadiusOuterRing() {return fOuterRadius;} + Float_t GetShiftX() { return fShiftX;} + Float_t GetShiftY() { return fShiftY;} + Float_t GetDetectorWhereX() { return fXcoord;} + Float_t GetDetectorWhereY() { return fYcoord;} + Float_t GetIntersectionX() {return fIntersectionX;} + Float_t GetIntersectionY() {return fIntersectionY;} + Float_t GetThetaOfRing() { return fThetaOfRing;} + Float_t GetAreaOfRing() {return fAreaOfRing;} + Float_t GetPortionOfRing() {return fPortionOfRing;} + Float_t GetHoughArea() { return fHoughArea;} + Float_t GetPhotonEta() { return fPhotonEta[fPhotonIndex];} + Float_t GetPhotonWeight() { return fPhotonWeight[fPhotonIndex];} + Float_t GetHoughRMS() { return fHoughRMS;} + Float_t* GetCandidatePhotonX() { return fCandidatePhotonX;} + Float_t* GetCandidatePhotonY() { return fCandidatePhotonY;} + Float_t GetHoughPhotonsNorm() { return fHoughPhotonsNorm;} + Float_t GetFittedTrackTheta() { return fFittedTrackTheta;} + Float_t GetFittedTrackPhi() { return fFittedTrackPhi;} + Float_t GetFittedThetaCerenkov() { return fFittedThetaCerenkov;} + Float_t GetEstimationOfTheta() { return fEstimationOfTheta;} + Float_t GetEstimationOfThetaRMS() { return fEstimationOfThetaRMS;} + void SetEventVertexZ(Float_t EventVertZ) { fEventVertZ = EventVertZ;} + void SetEventMultiplicity(Float_t EventMultiplicity) { fEventMultiplicity = EventMultiplicity;} + void SetPhotonEnergy(Float_t PhotonEnergy) { fPhotonEnergy = PhotonEnergy;} + void SetFreonRefractiveIndex() {fFreonRefractiveIndex = fFreonScaleFactor*(1.177+0.0172*fPhotonEnergy);} + void SetQuartzRefractiveIndex() {fQuartzRefractiveIndex = sqrt(1+(46.411/(113.763556-TMath::Power(fPhotonEnergy,2)))+(228.71/(328.51563-TMath::Power(fPhotonEnergy,2))));} + void SetGasRefractiveIndex() { fGasRefractiveIndex = 1.;} + void SetFreonScaleFactor(Float_t FreonScaleFactor) {fFreonScaleFactor = FreonScaleFactor;} + void SetEmissionPoint(Float_t LengthEmissionPoint) { fLengthEmissionPoint = LengthEmissionPoint;} + void SetMassHypotesis(Float_t mass) {fMassHypotesis = mass;} + void SetBetaOfParticle() { fTrackBeta = fTrackMomentum/sqrt(TMath::Power(fTrackMomentum,2)+TMath::Power(fMassHypotesis,2));} + void SetEntranceX(Float_t Xtoentr) { fXtoentr = Xtoentr;} + void SetEntranceY(Float_t Ytoentr) { fYtoentr = Ytoentr;} + void SetThetaPhotonInTRS(Float_t Theta) {fThetaPhotonInTRS = Theta;} + void SetPhiPhotonInTRS(Float_t Phi) {fPhiPhotonInTRS = Phi;} + void SetThetaPhotonInDRS(Float_t Theta) {fThetaPhotonInDRS = Theta;} + void SetPhiPhotonInDRS(Float_t Phi) {fPhiPhotonInDRS = Phi;} + void SetThetaAtQuartz(Float_t ThetaAtQuartz) {fThetaAtQuartz = ThetaAtQuartz;} + void SetPhiPoint(Float_t PhiPoint){ fPhiPoint = PhiPoint;} + void SetXCoordOfEmission(Float_t XEmiss) {fXEmiss = XEmiss;} + void SetYCoordOfEmission(Float_t YEmiss) {fYEmiss = YEmiss;} + void SetXPointOnCathode(Float_t PhotonLimitX) { fPhotonLimitX = PhotonLimitX;} + void SetYPointOnCathode(Float_t PhotonLimitY) { fPhotonLimitY = PhotonLimitY;} + void SetXInnerRing(Float_t XInner) {fXInner = XInner;} + void SetYInnerRing(Float_t YInner) {fYInner = YInner;} + void SetRadiusInnerRing(Float_t InnerRadius) {fInnerRadius = InnerRadius;} + void SetXOuterRing(Float_t XOuter) {fXOuter = XOuter;} + void SetYOuterRing(Float_t YOuter) {fYOuter = YOuter;} + void SetRadiusOuterRing(Float_t OuterRadius) {fOuterRadius = OuterRadius;} + void SetThetaCerenkov(Float_t ThetaCer) {fThetaCerenkov = ThetaCer;} + void SetThetaPhotonCerenkov(Float_t ThetaPhotCer) {fThetaPhotonCerenkov = ThetaPhotCer;} + void SetTrackMomentum(Float_t TrackMomentum) {fTrackMomentum = TrackMomentum;} + void SetTrackEta(Float_t TrackEta) {fTrackEta = TrackEta;} + void SetTrackTheta(Float_t TrackTheta) { fTrackTheta = TrackTheta;} + void SetTrackPhi(Float_t TrackPhi) { fTrackPhi = TrackPhi;} + void SetTrackPt(Float_t TrackPt) { fTrackPt = TrackPt;} + void SetTrackCharge(Int_t TrackCharge) { fTrackCharge = TrackCharge;} + void SetTrackTPCLastZ(Float_t TrackTPCLastZ) { fTrackTPCLastZ = TrackTPCLastZ;} + void SetMinDist(Float_t MinDist) { fMinDist = MinDist;} + void SetShiftX(Float_t ShiftX) { fShiftX = ShiftX;} + void SetShiftY(Float_t ShiftY) { fShiftY = ShiftY;} + void SetDetectorWhereX(Float_t Xcoord) { fXcoord = Xcoord;} + void SetDetectorWhereY(Float_t Ycoord) { fYcoord = Ycoord;} + void SetIntersectionX(Float_t IntersectionX) { fIntersectionX = IntersectionX;} + void SetIntersectionY(Float_t IntersectionY) { fIntersectionY = IntersectionY;} + void SetThetaOfRing(Float_t ThetaOfRing) { fThetaOfRing = ThetaOfRing;} + void SetAreaOfRing(Float_t AreaOfRing) { fAreaOfRing = AreaOfRing;} + void SetPortionOfRing(Float_t PortionOfRing) { fPortionOfRing = PortionOfRing;} + void SetHoughArea(Float_t HoughArea) { fHoughArea = HoughArea;} + void SetPhotonsNumber(Int_t PhotonsNumber) { fPhotonsNumber = PhotonsNumber;} + void SetPhotonIndex(Int_t PhotonIndex) { fPhotonIndex = PhotonIndex;} + void SetPhotonEta(Float_t PhotonEta) { fPhotonEta[fPhotonIndex] = PhotonEta;} + void SetPhotonFlag(Int_t PhotonFlag) { fPhotonFlag[fPhotonIndex] = PhotonFlag;} + void SetPhotonWeight(Float_t PhotonWeight) { fPhotonWeight[fPhotonIndex] = PhotonWeight;} + void SetHoughRMS(Float_t HoughRMS) { fHoughRMS = HoughRMS;} + void SetMipIndex(Int_t MipIndex) { fMipIndex = MipIndex;} + void SetTrackIndex(Int_t TrackIndex) { fTrackIndex = TrackIndex;} + void SetCandidatePhotonX(Float_t *CandidatePhotonX) { fCandidatePhotonX = CandidatePhotonX;} + void SetCandidatePhotonY(Float_t *CandidatePhotonY) { fCandidatePhotonY = CandidatePhotonY;} + void SetCandidatePhotonsNumber(Int_t CandidatePhotonsNumber) { fCandidatePhotonsNumber = CandidatePhotonsNumber;} + void SetHoughPhotons(Int_t HoughPhotons) { fHoughPhotons = HoughPhotons;} + void SetHoughPhotonsNorm(Float_t HoughPhotonsNorm) { fHoughPhotonsNorm = HoughPhotonsNorm;} + void SetFittedTrackTheta(Float_t FittedTrackTheta) { fFittedTrackTheta = FittedTrackTheta;} + void SetFittedTrackPhi(Float_t FittedTrackPhi) { fFittedTrackPhi = FittedTrackPhi;} + void SetFittedThetaCerenkov(Float_t FittedThetaCerenkov) { fFittedThetaCerenkov = FittedThetaCerenkov;} + void SetFittedHoughPhotons(Int_t FittedHoughPhotons) { fFittedHoughPhotons = FittedHoughPhotons;} + void SetEstimationOfTheta(Float_t EstimationOfTheta) { fEstimationOfTheta = EstimationOfTheta;} + void SetEstimationOfThetaRMS(Float_t EstimationOfThetaRMS) { fEstimationOfThetaRMS = EstimationOfThetaRMS;} + void FindBetaFromTheta(Float_t ThetaCerenkov) {fTrackBeta = 1/(fFreonRefractiveIndex*cos(ThetaCerenkov));} + Float_t SnellAngle(Float_t n1, Float_t n2, Float_t theta1); + Float_t FromEmissionToCathode(); +public: Int_t fNumEtaPhotons; // Number of photons Int_t fEtaFlag[3000]; // flag for good photons Float_t fEtaPhotons[3000]; // Cerenkov angle each photon @@ -344,10 +188,88 @@ class AliRICHRecon : public TTask { Float_t fThetaCerenkov; // Theta angle for Hough Float_t fWeightThetaCerenkov; // Theta Cerenkov angle weighted Float_t fThetaPeakPos; // Peak position - - -ClassDef(AliRICHRecon,0) - +protected: + AliRICH* fRich; + Int_t fTrackCharge; + Int_t fMipIndex; + Int_t fTrackIndex; + Int_t fPhotonsNumber; + Int_t fPhotonIndex; + Int_t fPhotonFlag[3000]; + Int_t fCandidatePhotonsNumber; + Int_t fHoughPhotons; + Int_t fFittedHoughPhotons; + Float_t fEventVertZ; + Float_t fEventMultiplicity; + Float_t fTrackTheta; + Float_t fTrackPhi; + Float_t fTrackMomentum; + Float_t fTrackEta; + Float_t fTrackPt; + Float_t fTrackTPCLastZ; + Float_t fMinDist; + Float_t fTrackBeta; + Float_t fXtoentr; + Float_t fYtoentr; + Float_t fThetaPhotonInTRS; + Float_t fPhiPhotonInTRS; + Float_t fThetaPhotonInDRS; + Float_t fPhiPhotonInDRS; + Float_t fThetaAtQuartz; + Float_t fPhiPoint; + Float_t fXEmiss; + Float_t fYEmiss; + Float_t fXInner; + Float_t fYInner; + Float_t fXOuter; + Float_t fYOuter; + Float_t fInnerRadius; + Float_t fOuterRadius; + Float_t fPhotonEnergy; + Float_t fFreonRefractiveIndex; + Float_t fQuartzRefractiveIndex; + Float_t fGasRefractiveIndex; + Float_t fFreonScaleFactor; + Float_t fLengthEmissionPoint; + Float_t fPhotonLimitX; + Float_t fPhotonLimitY; + Float_t fDistanceFromCluster; + Float_t fMassHypotesis; + Float_t fCerenkovAnglePad; + Float_t fThetaPhotonCerenkov; + Float_t fShiftX; + Float_t fShiftY; + Float_t fXcoord; + Float_t fYcoord; + Float_t fIntersectionX; + Float_t fIntersectionY; + Float_t fThetaOfRing; + Float_t fAreaOfRing; + Float_t fPortionOfRing; + Float_t fHoughArea; + Float_t fPhotonEta[3000]; + Float_t fPhotonWeight[3000]; + Float_t fHoughRMS; + Float_t* fCandidatePhotonX; + Float_t* fCandidatePhotonY; + Float_t fHoughPhotonsNorm; + Float_t fFittedTrackTheta; + Float_t fFittedTrackPhi; + Float_t fFittedThetaCerenkov; + Float_t fEstimationOfTheta; + Float_t fEstimationOfThetaRMS; + Int_t fThetaBin; Float_t fThetaMin,fThetaMax; + Float_t fXmin,fXmax,fYmin,fYmax; + TFile *fOutFile; + TNtuple *fNtuple; + TCanvas *fDisplay; + Int_t fNrings; //current number of reconstructed rings + Bool_t fDebug; + Bool_t fIsDISPLAY; + Bool_t fIsWEIGHT; + Bool_t fIsBACKGROUND; + Bool_t fIsMINIMIZER; + ClassDef(AliRICHRecon,0) }; #endif // #ifdef AliRICHRecon_cxx -- 2.43.0