X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RICH%2FAliRICH.cxx;h=924506f3031faeaa11460c12ad10567ab2b61cb6;hb=f701e8863eb28607018eb7ae3a8f151bee52798b;hp=5c39cf4f6446fd8bb8b2cc619f011539b2b87c9b;hpb=a3d71079810d242d0cb64a6c84a3512c19b6c1b2;p=u%2Fmrichter%2FAliRoot.git diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index 5c39cf4f644..924506f3031 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -15,6 +15,48 @@ /* $Log$ + + Revision 1.60 2002/10/22 16:28:21 alibrary + Introducing Riostream.h + + Revision 1.59 2002/10/14 14:57:31 hristov + Merging the VirtualMC branch to the main development branch (HEAD) + + Revision 1.58.6.1 2002/06/10 15:12:46 hristov + Merged with v3-08-02 + + Revision 1.58 2001/11/14 09:49:37 dibari + Use debug methods + + Revision 1.57 2001/11/09 17:29:31 dibari + Setters fro models moved to header + + Revision 1.56 2001/11/02 15:37:25 hristov + Digitizer class created. Code cleaning and bug fixes (J.Chudoba) + + Revision 1.55 2001/10/23 13:03:35 hristov + The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba) + + Revision 1.54 2001/09/07 08:38:10 hristov + Pointers initialised to 0 in the default constructors + + Revision 1.53 2001/08/30 09:51:23 hristov + The operator[] is replaced by At() or AddAt() in case of TObjArray. + + Revision 1.52 2001/05/16 14:57:20 alibrary + New files for folders and Stack + + Revision 1.51 2001/05/14 10:18:55 hristov + Default arguments declared once + + Revision 1.50 2001/05/10 14:44:16 jbarbosa + Corrected some overlaps (thanks I. Hrivnacovna). + + Revision 1.49 2001/05/10 12:23:49 jbarbosa + Repositioned the RICH modules. + Eliminated magic numbers. + Incorporated diagnostics (from macros). + Revision 1.48 2001/03/15 10:35:00 jbarbosa Corrected bug in MakeBranch (was using a different version of STEER) @@ -140,10 +182,10 @@ #include #include #include -//#include +#include #include -#include +#include #include #include "AliRICH.h" @@ -160,25 +202,35 @@ #include "AliRICHHitMapA1.h" #include "AliRICHClusterFinder.h" #include "AliRICHMerger.h" +#include "AliRICHDigitizer.h" #include "AliRun.h" +#include "AliRunDigitizer.h" #include "AliMC.h" #include "AliMagF.h" #include "AliConst.h" #include "AliPDG.h" #include "AliPoints.h" -#include "AliCallf77.h" -// Static variables for the pad-hit iterator routines -static Int_t sMaxIterPad=0; + + +static Int_t sMaxIterPad=0; // Static variables for the pad-hit iterator routines static Int_t sCurIterPad=0; ClassImp(AliRICH) //___________________________________________ +// RICH manager class +//Begin_Html +/* + +*/ +//End_Html + AliRICH::AliRICH() { -// Default constructor for RICH manager class +// Default ctor should not contain any new operators + cout<\n"; // no way to control it as ctor is called before call to SetDebugXXXX() fIshunt = 0; fHits = 0; @@ -191,33 +243,28 @@ AliRICH::AliRICH() fRawClusters = 0; fChambers = 0; fCerenkovs = 0; - for (Int_t i=0; i<7; i++) - { + for (Int_t i=0; i<7; i++){ fNdch[i] = 0; fNrawch[i] = 0; fNrechits1D[i] = 0; fNrechits3D[i] = 0; - } + } fFileName = 0; -} + fMerger = 0; +}//AliRICH::AliRICH() -//___________________________________________ AliRICH::AliRICH(const char *name, const char *title) : AliDetector(name,title) { -//Begin_Html -/* - -*/ -//End_Html - +// Named ctor + cout<\n"; // no way to control it as ctor is called before call to SetDebugXXXX() + fHits = new TClonesArray("AliRICHHit",1000 ); gAlice->AddHitList(fHits); fSDigits = new TClonesArray("AliRICHSDigit",100000); fCerenkovs = new TClonesArray("AliRICHCerenkov",1000); gAlice->AddHitList(fCerenkovs); - //gAlice->AddHitList(fHits); fNSDigits = 0; fNcerenkovs = 0; fIshunt = 0; @@ -232,7 +279,8 @@ AliRICH::AliRICH(const char *name, const char *title) Int_t i; for (i=0; iAddAt(new TClonesArray("AliRICHDigit",10000), i); fNdch[i]=0; } @@ -242,17 +290,20 @@ AliRICH::AliRICH(const char *name, const char *title) //printf("Created fRwClusters with adress:%p",fRawClusters); for (i=0; iAddAt(new TClonesArray("AliRICHRawCluster",10000), i); fNrawch[i]=0; } //fNrechits = new Int_t[kNCH]; for (i=0; iAddAt(new TClonesArray("AliRICHRecHit1D",1000), i); } for (i=0; iAddAt(new TClonesArray("AliRICHRecHit3D",1000), i); } //printf("Created fRecHits with adress:%p",fRecHits); @@ -264,19 +315,19 @@ AliRICH::AliRICH(const char *name, const char *title) (*fChambers)[i] = new AliRICHChamber();*/ fFileName = 0; + fMerger = 0; } AliRICH::AliRICH(const AliRICH& RICH) { -// Copy Constructor +// Copy ctor } -//___________________________________________ AliRICH::~AliRICH() { - -// Destructor of RICH manager class +// Dtor of RICH manager class + if(IsDebugStart()) cout<\n"; fIshunt = 0; delete fHits; @@ -311,20 +362,20 @@ AliRICH::~AliRICH() //_____________________________________________________________________________ Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res) { -// // Calls the charge disintegration method of the current chamber and adds -// the simulated cluster to the root treee -// - Int_t clhits[5]; - Float_t newclust[4][500]; - Int_t nnew; +// the simulated cluster to the root tree + if(IsDebugHit()||IsDebugDigit()) cout<\n"; + + Int_t clhits[5]; + Float_t newclust[4][500]; + Int_t nnew; // // Integrated pulse height on chamber clhits[0]=fNhits+1; - ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res); + ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, nnew, newclust, res); Int_t ic=0; // @@ -347,21 +398,21 @@ Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol } } - if (gAlice->TreeS()) - { + if (gAlice->TreeS()){ gAlice->TreeS()->Fill(); gAlice->TreeS()->Write(0,TObject::kOverwrite); //printf("Filled SDigits...\n"); - } + } -return nnew; -} -//___________________________________________ + return nnew; +}//Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res) + void AliRICH::Hits2SDigits() { - // Dummy: sdigits are created during transport. -// Called from alirun. +// Called from alirun. + if(IsDebugHit()||IsDebugDigit()) cout<\n"; + int nparticles = gAlice->GetNtrack(); cout << "Particles (RICH):" <\n"; - AliRICHChamber* iChamber; + //AliRICHChamber* iChamber; - printf("Generating tresholds...\n"); + //printf("Generating tresholds...\n"); - for(Int_t i=0;i<7;i++) { - iChamber = &(Chamber(i)); - iChamber->GenerateTresholds(); - } + //for(Int_t i=0;i<7;i++) { + //iChamber = &(Chamber(i)); + //iChamber->GenerateTresholds(); + //} - int nparticles = gAlice->GetNtrack(); - if (nparticles > 0) - { - if (fMerger) { - fMerger->Init(); - fMerger->Digitise(nev,flag); - } - } - //Digitise(nev,flag); + //int nparticles = gAlice->GetNtrack(); + //cout << "Particles (RICH):" <Init(); + //fMerger->Digitise(nev,flag); + + AliRunDigitizer * manager = new AliRunDigitizer(1,1); + manager->SetInputStream(0,"galice.root"); + //AliRICHDigitizer *dRICH = new AliRICHDigitizer(manager); + manager->Exec("deb"); + } //___________________________________________ void AliRICH::SDigits2Digits() { - -// -// Generate digits -// Called from alirun, single event only. - - AliRICHChamber* iChamber; - - printf("Generating tresholds...\n"); - - for(Int_t i=0;i<7;i++) { - iChamber = &(Chamber(i)); - iChamber->GenerateTresholds(); - } - - int nparticles = gAlice->GetNtrack(); - cout << "Particles (RICH):" < 0) - { - if (fMerger) { - fMerger->Init(); - fMerger->Digitise(0,0); - } - } + SDigits2Digits(0,0); } //___________________________________________ void AliRICH::Digits2Reco() { - // Generate clusters -// Called from alirun, single event only. +// Called from alirun, single event only. + if(IsDebugDigit()||IsDebugReco()) cout<\n"; + int nparticles = gAlice->GetNtrack(); cout << "Particles (RICH):" <\n"; TClonesArray &lhits = *fHits; new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits); } -//_____________________________________________________________________________ + void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs) { - -// -// Adds a RICH cerenkov hit to the Cerenkov Hits list -// +// Adds a RICH cerenkov hit to the Cerenkov Hits list + if(IsDebugHit()) cout<\n"; TClonesArray &lcerenkovs = *fCerenkovs; new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs); - //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs); } -//___________________________________________ -void AliRICH::AddSDigit(Int_t *clhits) -{ -// -// Add a RICH pad hit to the list -// +void AliRICH::AddSDigit(Int_t *aSDigit) +{ +// Adds the current S digit to the RICH list of S digits + if(IsDebugDigit()) cout<\n"; - //printf("fsdigits:%p, data: %d\n",fSDigits,clhits[2]); TClonesArray &lSDigits = *fSDigits; - new(lSDigits[fNSDigits++]) AliRICHSDigit(clhits); + new(lSDigits[fNSDigits++]) AliRICHSDigit(aSDigit); } -//_____________________________________________________________________________ + + void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits) { +// Add a RICH digit to the list + if(IsDebugDigit()) cout<\n"; - // - // Add a RICH digit to the list - // - - //printf("fdigits:%p, data: %d\n",((TClonesArray*)(*fDchambers)[id]),digits[0]); - TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]); - new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits); + TClonesArray &ldigits = *((TClonesArray*)fDchambers->At(id)); + new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits); } -//_____________________________________________________________________________ void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c) { - // - // Add a RICH digit to the list - // +// Add a RICH digit to the list + + if(IsDebugStart()) + cout<\n"; - TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]); + //PH TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]); + TClonesArray &lrawcl = *((TClonesArray*)fRawClusters->At(id)); new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c); } @@ -502,20 +529,18 @@ void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *pa // Add a RICH reconstructed hit to the list // - TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]); + //PH TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]); + TClonesArray &lrec1D = *((TClonesArray*)fRecHits1D->At(id)); new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy); } //_____________________________________________________________________________ -void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit) +void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit, Float_t omega, Float_t theta, Float_t phi) { - - // - // Add a RICH reconstructed hit to the list - // +// Add a RICH reconstructed hit to the list - TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]); - new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit); + TClonesArray &lrec3D = *((TClonesArray*)fRecHits3D->At(id)); + new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit,omega,theta,phi); } //___________________________________________ @@ -538,7 +563,7 @@ void AliRICH::BuildGeometry() AliRICHGeometry* geometry; iChamber = &(pRICH->Chamber(0)); - segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0); + segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(); geometry=iChamber->GetGeometryModel(); new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15); @@ -802,7 +827,7 @@ void AliRICH::CreateGeometry() AliRICHChamber* iChamber; iChamber = &(pRICH->Chamber(0)); - segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0); + segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(); geometry=iChamber->GetGeometryModel(); Float_t distance; @@ -1106,10 +1131,10 @@ void AliRICH::CreateGeometry() // --- Places the detectors defined with GSVOLU // Place material inside RICH gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY"); - gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY"); - gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY"); - gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, -68.35 - 1.25, 0, "ONLY"); - gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 68.35 + 1.25, 0, "ONLY"); + gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY"); + gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY"); + gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY"); + gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 68.35 + 1.25, 0, "ONLY"); gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY"); @@ -1213,8 +1238,8 @@ void AliRICH::CreateGeometry() // PCB placing - gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 1.2, 0, "ONLY"); - gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 1.2, 0, "ONLY"); + gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 2.5, 0, "ONLY"); + gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY"); @@ -1813,7 +1838,7 @@ Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t ) } //___________________________________________ -void AliRICH::MakeBranch(Option_t* option, char *file) +void AliRICH::MakeBranch(Option_t* option, const char *file) { // Create Tree branches for the RICH. @@ -1831,14 +1856,14 @@ void AliRICH::MakeBranch(Option_t* option, char *file) if (cH) { sprintf(branchname,"%sCerenkov",GetName()); if (fCerenkovs && gAlice->TreeH()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ; + //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); } sprintf(branchname,"%sSDigits",GetName()); if (fSDigits && gAlice->TreeH()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ; + //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); //printf("Making branch %sSDigits in TreeH\n",GetName()); } @@ -1847,8 +1872,8 @@ void AliRICH::MakeBranch(Option_t* option, char *file) if (cS) { sprintf(branchname,"%sSDigits",GetName()); if (fSDigits && gAlice->TreeS()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ; + //TBranch* branch = MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); //printf("Making branch %sSDigits in TreeS\n",GetName()); } @@ -1863,8 +1888,8 @@ void AliRICH::MakeBranch(Option_t* option, char *file) for (i=0; iTreeD()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ; + //TBranch* branch = MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); //printf("Making Branch %sDigits%d\n",GetName(),i+1); } @@ -1883,8 +1908,8 @@ void AliRICH::MakeBranch(Option_t* option, char *file) for (i=0; iTreeR()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ; + //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); } } @@ -1894,16 +1919,15 @@ void AliRICH::MakeBranch(Option_t* option, char *file) for (i=0; iTreeR()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ; + //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); } } for (i=0; iTreeR()) { - //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ; - gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ; //branch->SetAutoDelete(kFALSE); } } @@ -2007,7 +2031,8 @@ void AliRICH::ResetDigits() // Reset number of digits and the digits array for this detector // for ( int i=0;iClear(); + //PH if (fDchambers && (*fDchambers)[i]) (*fDchambers)[i]->Clear(); + if (fDchambers && fDchambers->At(i)) fDchambers->At(i)->Clear(); if (fNdch) fNdch[i]=0; } } @@ -2019,7 +2044,8 @@ void AliRICH::ResetRawClusters() // Reset number of raw clusters and the raw clust array for this detector // for ( int i=0;iClear(); + //PH if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear(); + if (fRawClusters->At(i)) ((TClonesArray*)fRawClusters->At(i))->Clear(); if (fNrawch) fNrawch[i]=0; } } @@ -2032,7 +2058,8 @@ void AliRICH::ResetRecHits1D() // for ( int i=0;iClear(); + //PH if ((*fRecHits1D)[i]) ((TClonesArray*)(*fRecHits1D)[i])->Clear(); + if (fRecHits1D->At(i)) ((TClonesArray*)fRecHits1D->At(i))->Clear(); if (fNrechits1D) fNrechits1D[i]=0; } } @@ -2045,59 +2072,16 @@ void AliRICH::ResetRecHits3D() // for ( int i=0;iClear(); + //PH if ((*fRecHits3D)[i]) ((TClonesArray*)(*fRecHits3D)[i])->Clear(); + if (fRecHits3D->At(i)) ((TClonesArray*)fRecHits3D->At(i))->Clear(); if (fNrechits3D) fNrechits3D[i]=0; } } -//___________________________________________ -void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry) -{ - -// -// Setter for the RICH geometry model -// - - - ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry); -} - -//___________________________________________ -void AliRICH::SetSegmentationModel(Int_t id, AliSegmentation *segmentation) -{ - -// -// Setter for the RICH segmentation model -// - - ((AliRICHChamber*) (*fChambers)[id])->SetSegmentationModel(segmentation); -} - -//___________________________________________ -void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response) -{ - -// -// Setter for the RICH response model -// - - ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response); -} - -void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst) -{ - -// -// Setter for the RICH reconstruction model (clusters) -// - - ((AliRICHChamber*) (*fChambers)[id])->SetReconstructionModel(reconst); -} //___________________________________________ void AliRICH::StepManager() { - // Full Step Manager Int_t copy, id; @@ -2115,7 +2099,7 @@ void AliRICH::StepManager() Float_t localTheta,localPhi; Float_t theta,phi; Float_t destep, step; - Float_t ranf[2]; + Double_t ranf[2]; Int_t nPads; Float_t coscerenkov; static Float_t eloss, xhit, yhit, tlength; @@ -2214,7 +2198,8 @@ void AliRICH::StepManager() Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1])); Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi); - gMC->Rndm(ranf, 1); + //gMC->Rndm(ranf, 1); + gMC->GetRandom()->RndmArray(1,ranf); //printf("grid calculation:%f\n",t); if (ranf[0] > t) { gMC->StopTrack(); @@ -2242,7 +2227,8 @@ void AliRICH::StepManager() Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1])); Float_t t = Fresnel(ckovEnergy*1e9,cophi,1); - gMC->Rndm(ranf, 1); + //gMC->Rndm(ranf, 1); + gMC->GetRandom()->RndmArray(1,ranf); if (ranf[0] < t) { gMC->StopTrack(); ckovData[13] = 6; @@ -2344,8 +2330,19 @@ void AliRICH::StepManager() Double_t rt = TMath::Sqrt(tc); theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg; phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg; - gMC->Gmtod(pos,localPos,1); - gMC->Gmtod(mom,localMom,2); + + gMC->CurrentVolOffID(2,copy); + vol[0]=copy; + idvol=vol[0]-1; + + + //gMC->Gmtod(pos,localPos,1); + + Chamber(idvol).GlobaltoLocal(pos,localPos); + + //gMC->Gmtod(mom,localMom,2); + + Chamber(idvol).GlobaltoLocal(mom,localMom); gMC->CurrentVolOffID(2,copy); vol[0]=copy; @@ -2360,7 +2357,8 @@ void AliRICH::StepManager() printf("Feedbacks:%d\n",fFeedbacks); }*/ - ((AliRICHChamber*) (*fChambers)[idvol]) + //PH ((AliRICHChamber*) (*fChambers)[idvol]) + ((AliRICHChamber*)fChambers->At(idvol)) ->SigGenInit(localPos[0], localPos[2], localPos[1]); if(idvolTrackPid(); // particle type @@ -2400,9 +2398,9 @@ void AliRICH::StepManager() mom[0] = current->Px(); mom[1] = current->Py(); mom[2] = current->Pz(); - Float_t mipPx = mipHit->fMomX; - Float_t mipPy = mipHit->fMomY; - Float_t mipPz = mipHit->fMomZ; + Float_t mipPx = mipHit->MomX(); + Float_t mipPy = mipHit->MomY(); + Float_t mipPz = mipHit->MomZ(); Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2]; Float_t rt = TMath::Sqrt(r); @@ -2476,8 +2474,14 @@ void AliRICH::StepManager() mom[1]=momentum(1); mom[2]=momentum(2); mom[3]=momentum(3); - gMC->Gmtod(pos,localPos,1); - gMC->Gmtod(mom,localMom,2); + + //gMC->Gmtod(pos,localPos,1); + + Chamber(idvol).GlobaltoLocal(pos,localPos); + + //gMC->Gmtod(mom,localMom,2); + + Chamber(idvol).GlobaltoLocal(mom,localMom); ipart = gMC->TrackPid(); // @@ -2528,7 +2532,8 @@ void AliRICH::StepManager() if(idvolAt(idvol)) ->SigGenInit(localPos[0], localPos[2], localPos[1]); } } @@ -2570,10 +2575,12 @@ void AliRICH::StepManager() // defined by the segmentation // model (boundary crossing conditions) } else if - (((AliRICHChamber*) (*fChambers)[idvol]) + //PH (((AliRICHChamber*) (*fChambers)[idvol]) + (((AliRICHChamber*)fChambers->At(idvol)) ->SigGenCond(localPos[0], localPos[2], localPos[1])) { - ((AliRICHChamber*) (*fChambers)[idvol]) + //PH ((AliRICHChamber*) (*fChambers)[idvol]) + ((AliRICHChamber*)fChambers->At(idvol)) ->SigGenInit(localPos[0], localPos[2], localPos[1]); if (eloss > 0) { @@ -2597,7 +2604,7 @@ void AliRICH::StepManager() } /*************************************************End of MIP treatment**************************************/ //} -} +}//void AliRICH::StepManager() void AliRICH::FindClusters(Int_t nev,Int_t lastEntry) { @@ -2609,7 +2616,8 @@ void AliRICH::FindClusters(Int_t nev,Int_t lastEntry) gAlice->ResetDigits(); gAlice->TreeD()->GetEvent(0); for (Int_t ich=0;ichAt(ich); TClonesArray *pRICHdigits = this->DigitsAddress(ich); if (pRICHdigits == 0) continue; @@ -2659,9 +2667,9 @@ AliRICHSDigit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters ) // Return the address of the first sdigit for hit TClonesArray *theClusters = clusters; Int_t nclust = theClusters->GetEntriesFast(); - if (nclust && hit->fPHlast > 0) { - sMaxIterPad=Int_t(hit->fPHlast); - sCurIterPad=Int_t(hit->fPHfirst); + if (nclust && hit->PHlast() > 0) { + sMaxIterPad=Int_t(hit->PHlast()); + sCurIterPad=Int_t(hit->PHfirst()); return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1); } else { return 0; @@ -2689,7 +2697,7 @@ AliRICH& AliRICH::operator=(const AliRICH& rhs) } -void AliRICH::DiagnosticsFE(Int_t evNumber1=0,Int_t evNumber2=0) +void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2) { Int_t NpadX = 162; // number of pads on X @@ -2778,12 +2786,12 @@ void AliRICH::DiagnosticsFE(Int_t evNumber1=0,Int_t evNumber2=0) //Float_t x = mHit->X(); // x-pos of hit //Float_t y = mHit->Z(); // y-pos //Float_t z = mHit->Y(); - //Float_t phi = mHit->fPhi; //Phi angle of incidence - Float_t theta = mHit->fTheta; //Theta angle of incidence + //Float_t phi = mHit->Phi(); //Phi angle of incidence + Float_t theta = mHit->Theta(); //Theta angle of incidence Float_t px = mHit->MomX(); Float_t py = mHit->MomY(); Int_t index = mHit->Track(); - Int_t particle = (Int_t)(mHit->fParticle); + Int_t particle = (Int_t)(mHit->Particle()); Float_t R; Float_t PTfinal; Float_t PTvertex; @@ -3001,6 +3009,10 @@ void AliRICH::DiagnosticsFE(Int_t evNumber1=0,Int_t evNumber2=0) } // } + + TStyle *mystyle=new TStyle("Plain","mystyle"); + mystyle->SetPalette(1,0); + mystyle->cd(); //Create canvases, set the view range, show histograms @@ -3207,7 +3219,7 @@ void AliRICH::DiagnosticsFE(Int_t evNumber1=0,Int_t evNumber2=0) //_________________________________________________________________________________________________ -void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1=0,Int_t evNumber2=0) +void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2) { AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); @@ -3215,7 +3227,7 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); AliRICHChamber* chamber; chamber = &(pRICH->Chamber(0)); - segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel(0); + segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel(); Int_t NpadX = segmentation->Npx(); // number of pads on X Int_t NpadY = segmentation->Npy(); // number of pads on Y @@ -3232,6 +3244,11 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); Int_t xmax= NpadX/2; Int_t ymin= -NpadY/2; Int_t ymax= NpadY/2; + + Float_t PTfinal = 0; + Int_t pionCount = 0; + Int_t kaonCount = 0; + Int_t protonCount = 0; TH2F *feedback = 0; TH2F *mip = 0; @@ -3240,17 +3257,17 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); TH1F *hitsX = 0; TH1F *hitsY = 0; - TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-30,30,150,-50,10); + TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-25,25,150,-45,5); if (diaglevel == 1) { printf("Single Ring Hits\n"); - feedback = new TH2F("feedback","Feedback hit distribution",150,-30,30,150,-50,10); - mip = new TH2F("mip","Mip hit distribution",150,-30,30,150,-50,10); - cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-30,30,150,-50,10); - h = new TH2F("h","Detector hit distribution",150,-30,30,150,-50,10); - hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-30,30); - hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,10); + feedback = new TH2F("feedback","Feedback hit distribution",150,-20,20,150,-35,5); + mip = new TH2F("mip","Mip hit distribution",150,-20,20,150,-35,5); + cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-20,20,150,-35,5); + h = new TH2F("h","Detector hit distribution",150,-20,20,150,-35,5); + hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-50,50); + hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,50); } else { @@ -3275,7 +3292,7 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax); TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.); - TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",200,.3,1); + TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",100,.35,.8); TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1); TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.); TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.); @@ -3295,15 +3312,23 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.); TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.); TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.); - TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,0,360); - TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence",100,0,15); - TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.5,1); - TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",200,0,15); - TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",200,0,360); - TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.3,1); - TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",200,.3,1); + TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",50,0,360); + TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of theta angle of incidence",50,0,15); + TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",50,.5,1); + TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",100,0,15); + TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",100,0,360); + TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",100,.35,.8); + TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",100,.35,.8); TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30); TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.); + TH2F *identification = new TH2F("identification","Particle Identification",100,1,5,100,0,.8); + TH1F *OriginalOmega = new TH1F("Original Omega","Cerenkov angle per track",100,.35,.8); + TH1F *OriginalPhi = new TH1F("Original Phi","Distribution of phi angle of incidence per track",100,0,360); + TH1F *OriginalTheta = new TH1F("Original Theta","Distribution of theta angle per track",100,0,15); + TH1F *OmegaError = new TH1F("Omega Error","Difference between original an reconstructed cerenkov angle",100,0,.2); + TH1F *PhiError = new TH1F("Phi Error","Difference between original an reconstructed phi angle",100,0,360); + TH1F *ThetaError = new TH1F("Theta Error","Difference between original an reconstructed phi angle",100,0,15); + // Start loop over events @@ -3318,6 +3343,17 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); Int_t feed=0; Int_t padmip=0; Float_t x=0,y=0; + + Float_t chiSquareOmega = 0; + Float_t chiSquareTheta = 0; + Float_t chiSquarePhi = 0; + + Float_t recEffEvent = 0; + Float_t recEffTotal = 0; + + Float_t trackglob[3]; + Float_t trackloc[3]; + for (Int_t i=0;i<100;i++) mothers[i]=0; @@ -3352,14 +3388,37 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); mHit; mHit=(AliRICHHit*)pRICH->NextHit()) { - //Int_t nch = mHit->fChamber; // chamber number - x = mHit->X(); // x-pos of hit - y = mHit->Z(); // y-pos - Float_t phi = mHit->fPhi; //Phi angle of incidence - Float_t theta = mHit->fTheta; //Theta angle of incidence + Int_t nch = mHit->Chamber(); // chamber number + trackglob[0] = mHit->X(); // x-pos of hit + trackglob[1] = mHit->Y(); + trackglob[2] = mHit->Z(); // y-pos of hit + //x = mHit->X(); // x-pos of hit + //y = mHit->Z(); // y-pos + Float_t phi = mHit->Phi(); //Phi angle of incidence + Float_t theta = mHit->Theta(); //Theta angle of incidence Int_t index = mHit->Track(); - Int_t particle = (Int_t)(mHit->fParticle); + Int_t particle = (Int_t)(mHit->Particle()); //Int_t freon = (Int_t)(mHit->fLoss); + Float_t px = mHit->MomX(); + Float_t py = mHit->MomY(); + + if (TMath::Abs(particle) < 10000000) + { + PTfinal=TMath::Sqrt(px*px + py*py); + //printf("PTfinal 0: %f\n",PTfinal); + } + + chamber = &(pRICH->Chamber(nch-1)); + + //printf("Nch:%d\n",nch); + + chamber->GlobaltoLocal(trackglob,trackloc); + + chamber->LocaltoGlobal(trackloc,trackglob); + + + x=trackloc[0]; + y=trackloc[2]; hitsX->Fill(x,(float) 1); hitsY->Fill(y,(float) 1); @@ -3413,7 +3472,7 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); h->Fill(x,y,(float) 1); //} //} - } + } Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast(); //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040) @@ -3424,16 +3483,33 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); totalphotonsevent->Fill(ncerenkovs,(float) 1); for (Int_t hit=0;hitCerenkovs()->UncheckedAt(hit); - //Int_t nchamber = cHit->fChamber; // chamber number + Int_t nchamber = cHit->fChamber; // chamber number Int_t index = cHit->Track(); //Int_t pindex = (Int_t)(cHit->fIndex); - Float_t cx = cHit->X(); // x-position - Float_t cy = cHit->Z(); // y-position + trackglob[0] = cHit->X(); // x-pos of hit + trackglob[1] = cHit->Y(); + trackglob[2] = cHit->Z(); // y-pos of hit + //Float_t cx = cHit->X(); // x-position + //Float_t cy = cHit->Z(); // y-position Int_t cmother = cHit->fCMother; // Index of mother particle Int_t closs = (Int_t)(cHit->fLoss); // How did the particle get lost? Float_t cherenkov = cHit->fCerenkovAngle; //production cerenkov angle - //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy); + chamber = &(pRICH->Chamber(nchamber-1)); + + //printf("Nch:%d\n",nch); + + chamber->GlobaltoLocal(trackglob,trackloc); + + chamber->LocaltoGlobal(trackloc,trackglob); + + + Float_t cx=trackloc[0]; + Float_t cy=trackloc[2]; + + //printf ("Cerenkov hit number %d/%d, X:%f, Y:%f\n",hit,ncerenkovs,cx,cy); + + //printf("Particle:%9d\n",index); TParticle *current = (TParticle*)gAlice->Particle(index); @@ -3496,8 +3572,8 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); Float_t mx = mipHit->X(); Float_t my = mipHit->Z(); //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my); - Float_t dx = cx - mx; - Float_t dy = cy - my; + Float_t dx = trackglob[0] - mx; + Float_t dy = trackglob[2] - my; //printf("Dx:%f, Dy:%f\n",dx,dy); Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy); //printf("Final radius:%f\n",final_radius); @@ -3588,28 +3664,115 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); if(nrechits3D) { - for (Int_t hit=0;hitRecHitsAddress3D(2)->UncheckedAt(hit); - Float_t r_omega = recHit3D->fOmega; // Cerenkov angle - Float_t r_theta = recHit3D->fTheta; // Theta angle of incidence - Float_t r_phi = recHit3D->fPhi; // Phi angle if incidence - Float_t meanradius = recHit3D->fMeanRadius; // Mean radius for reconstructed point - - //printf("rechit %f %f %f %f %f\n",recHit3D->fOmega,recHit3D->fTheta,recHit3D->fPhi, recHit3D->fX,recHit3D->fY); - - Omega3D->Fill(r_omega,(float) 1); - Theta->Fill(r_theta*180/TMath::Pi(),(float) 1); - Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1); - MeanRadius->Fill(meanradius,(float) 1); - } + recEffEvent = 0; + + //for (Int_t hit=0;hitRecHitsAddress3D(2)->UncheckedAt(track); + Float_t r_omega = recHit3D->fOmega; // Cerenkov angle + Float_t r_theta = recHit3D->fTheta; // Theta angle of incidence + Float_t r_phi = recHit3D->fPhi; // Phi angle if incidence + Float_t meanradius = recHit3D->fMeanRadius; // Mean radius for reconstructed point + Float_t originalOmega = recHit3D->fOriginalOmega; // Real Cerenkov angle + Float_t originalTheta = recHit3D->fOriginalTheta; // Real incidence angle + Float_t originalPhi = recHit3D->fOriginalPhi; // Real azimuthal angle + + + //correction to track cerenkov angle + originalOmega = (Float_t) ckovangle->GetMean(); + + if(diaglevel == 4) + { + printf("\nMean cerenkov angle: %f\n", originalOmega); + printf("Reconstructed cerenkov angle: %f\n",r_omega); + } + + Float_t omegaError = TMath::Abs(originalOmega - r_omega); + Float_t thetaError = TMath::Abs(originalTheta - r_theta); + Float_t phiError = TMath::Abs(originalPhi - r_phi); + + //chiSquareOmega += (omegaError/originalOmega)*(omegaError/originalOmega); + //chiSquareTheta += (thetaError/originalTheta)*(thetaError/originalTheta); + //chiSquarePhi += (phiError/originalPhi)*(phiError/originalPhi); + + if(TMath::Abs(omegaError) < 0.015) + recEffEvent += 1; + + + + //printf("rechit %f %f %f %f %f\n",recHit3D->fOmega,recHit3D->fTheta,recHit3D->fPhi, recHit3D->fX,recHit3D->fY); + + Omega3D->Fill(r_omega,(float) 1); + Theta->Fill(r_theta*180/TMath::Pi(),(float) 1); + Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1); + MeanRadius->Fill(meanradius,(float) 1); + identification->Fill(PTfinal, r_omega,1); + OriginalOmega->Fill(originalOmega, (float) 1); + OriginalTheta->Fill(originalTheta, (float) 1); + OriginalPhi->Fill(TMath::Abs(originalPhi), (float) 1); + OmegaError->Fill(omegaError, (float) 1); + ThetaError->Fill(thetaError, (float) 1); + PhiError->Fill(phiError, (float) 1); + + recEffEvent = recEffEvent; + recEffTotal += recEffEvent; + + Float_t pioncer = acos(sqrt((.139*.139+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285))); + Float_t kaoncer = acos(sqrt((.439*.439+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285))); + Float_t protoncer = acos(sqrt((.938*.938+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285))); + + Float_t piondist = TMath::Abs(r_omega - pioncer); + Float_t kaondist = TMath::Abs(r_omega - kaoncer); + Float_t protondist = TMath::Abs(r_omega - protoncer); + + if(diaglevel == 4) + { + if(pioncerr_omega) + { + if(kaondist>piondist) + { + printf("Identified as a PION!\n"); + pionCount += 1; + } + else + { + printf("Identified as a KAON!\n"); + kaonCount += 1; + } + } } + if(protoncerr_omega) + { + if(kaondist>protondist) + { + printf("Identified as a PROTON!\n"); + protonCount += 1; + } + else + { + printf("Identified as a KAON!\n"); + pionCount += 1; + } + } + if(protoncer>r_omega) + { + printf("Identified as a PROTON!\n"); + protonCount += 1; + } + + printf("\nReconstruction efficiency: %5.2f%%\n", recEffEvent*100); } } } + for (Int_t nmothers=0;nmothersFill(mothers[nmothers],(float) 1); - mother->Fill(mothers2[nmothers],(float) 1); - //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]); + totalphotonstrack->Fill(mothers[nmothers],(float) 1); + mother->Fill(mothers2[nmothers],(float) 1); + //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]); } clusev->Fill(nraw,(float) 1); @@ -3623,31 +3786,31 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); pads = 0; nraw=0; padmip=0; - - - + + + gAlice->ResetDigits(); //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries(); gAlice->TreeD()->GetEvent(0); - + if (diaglevel < 4) { - - + + TClonesArray *Digits = pRICH->DigitsAddress(2); Int_t ndigits = Digits->GetEntriesFast(); printf("Digits : %d\n",ndigits); padsev->Fill(ndigits,(float) 1); for (Int_t hit=0;hitUncheckedAt(hit); - Int_t qtot = dHit->fSignal; // charge - Int_t ipx = dHit->fPadX; // pad number on X - Int_t ipy = dHit->fPadY; // pad number on Y + Int_t qtot = dHit->Signal(); // charge + Int_t ipx = dHit->PadX(); // pad number on X + Int_t ipy = dHit->PadY(); // pad number on Y //printf("%d, %d\n",ipx,ipy); if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot); } } - + if (diaglevel == 5) { for (Int_t ich=0;ich<7;ich++) @@ -3659,11 +3822,11 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); if (ndigits) { for (Int_t hit=0;hitUncheckedAt(hit); - //Int_t nchamber = dHit->fChamber; // chamber number + //Int_t nchamber = dHit->GetChamber(); // chamber number //Int_t nhit = dHit->fHitNumber; // hit number - Int_t qtot = dHit->fSignal; // charge - Int_t ipx = dHit->fPadX; // pad number on X - Int_t ipy = dHit->fPadY; // pad number on Y + Int_t qtot = dHit->Signal(); // charge + Int_t ipx = dHit->PadX(); // pad number on X + Int_t ipy = dHit->PadY(); // pad number on Y //Int_t iqpad = dHit->fQpad; // charge per pad //Int_t rpad = dHit->fRSec; // R-position of pad //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy); @@ -3680,7 +3843,61 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); } } } + + if(diaglevel == 4) + { + + Stat_t omegaE; + Stat_t thetaE; + Stat_t phiE; + + Stat_t omegaO; + Stat_t thetaO; + Stat_t phiO; + + for(Int_t i=0;i<99;i++) + { + omegaE = OriginalOmega->GetBinContent(i); + if(omegaE != 0) + { + omegaO = Omega3D->GetBinContent(i); + chiSquareOmega += (TMath::Power(omegaE,2) - TMath::Power(omegaO,2))/omegaO; + } + + thetaE = OriginalTheta->GetBinContent(i); + if(thetaE != 0) + { + thetaO = Theta->GetBinContent(i); + chiSquareTheta += (TMath::Power(thetaE,2) - TMath::Power(thetaO,2))/thetaO; + } + + phiE = OriginalPhi->GetBinContent(i); + if(phiE != 0) + { + phiO = Phi->GetBinContent(i); + chiSquarePhi += (TMath::Power(phiE,2) - TMath::Power(phiO,2))/phiO; + } + + //printf(" o: %f t: %f p: %f\n", OriginalOmega->GetBinContent(i), OriginalTheta->GetBinContent(i),OriginalPhi->GetBinContent(i)); + + } + + + + printf("\nChi square test values: Omega - %f\n", chiSquareOmega); + printf(" Theta - %f\n", chiSquareTheta); + printf(" Phi - %f\n", chiSquarePhi); + printf("\nKolmogorov test values: Omega - %5.4f\n", Omega3D->KolmogorovTest(OriginalOmega)); + printf(" Theta - %5.4f\n", Theta->KolmogorovTest(OriginalTheta)); + printf(" Phi - %5.4f\n", Phi->KolmogorovTest(OriginalPhi)); + + recEffTotal = recEffTotal/evNumber2; + printf("\nTotal reconstruction efficiency: %5.2f%%\n", recEffTotal*100); + printf("\n Pions: %d\n Kaons: %d\n Protons:%d\n",pionCount, kaonCount, protonCount); + + } + //Create canvases, set the view range, show histograms @@ -3696,9 +3913,25 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); TCanvas *c10 = 0; TCanvas *c11 = 0; TCanvas *c12 = 0; + TCanvas *c13 = 0; + //TF1* expo = 0; //TF1* gaus = 0; + TStyle *mystyle=new TStyle("Plain","mystyle"); + mystyle->SetPalette(1,0); + //mystyle->SetTitleYSize(0.2); + //mystyle->SetStatW(0.19); + //mystyle->SetStatH(0.1); + //mystyle->SetStatFontSize(0.01); + //mystyle->SetTitleYSize(0.3); + mystyle->SetFuncColor(2); + //mystyle->SetOptStat(0111); + mystyle->SetDrawBorder(0); + mystyle->SetTitleBorderSize(0); + mystyle->SetOptFit(1111); + mystyle->cd(); + TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2); Int_t nrechits3D = RecHits3D->GetEntriesFast(); @@ -3711,7 +3944,7 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350); hc0->SetXTitle("ix (npads)"); - hc0->Draw("box"); + hc0->Draw("colz"); // c2 = new TCanvas("c2","Hits per type",100,100,600,700); @@ -3721,25 +3954,25 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); c2->cd(1); feedback->SetXTitle("x (cm)"); feedback->SetYTitle("y (cm)"); - feedback->Draw(); + feedback->Draw("colz"); c2->cd(2); //mip->SetFillColor(5); mip->SetXTitle("x (cm)"); mip->SetYTitle("y (cm)"); - mip->Draw(); + mip->Draw("colz"); c2->cd(3); //cerenkov->SetFillColor(5); cerenkov->SetXTitle("x (cm)"); cerenkov->SetYTitle("y (cm)"); - cerenkov->Draw(); + cerenkov->Draw("colz"); c2->cd(4); //h->SetFillColor(5); h->SetXTitle("x (cm)"); h->SetYTitle("y (cm)"); - h->Draw(); + h->Draw("colz"); c3 = new TCanvas("c3","Hits distribution",150,150,600,350); c3->Divide(2,1); @@ -3906,38 +4139,149 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); if(nrechits3D) { - c8 = new TCanvas("c8","3D reconstruction",50,50,1100,700); - c8->Divide(4,2); + c8 = new TCanvas("c8","3D reconstruction of Phi angle",50,50,300,1050); + c8->Divide(1,3); //c2->SetFillColor(42); + + // data per hit c8->cd(1); hitsPhi->SetFillColor(5); + if (evNumber2>10) + hitsPhi->Fit("gaus"); hitsPhi->Draw(); + + //data per track c8->cd(2); - hitsTheta->SetFillColor(5); - hitsTheta->Draw(); + OriginalPhi->SetFillColor(5); + if (evNumber2>10) + OriginalPhi->Fit("gaus"); + OriginalPhi->Draw(); + + //recontructed data c8->cd(3); - ckovangle->SetFillColor(5); - ckovangle->SetXTitle("angle (radians)"); - ckovangle->Draw(); - c8->cd(4); - radius->SetFillColor(5); - radius->SetXTitle("radius (cm)"); - radius->Draw(); - c8->cd(5); Phi->SetFillColor(5); + if (evNumber2>10) + Phi->Fit("gaus"); Phi->Draw(); - c8->cd(6); + + c9 = new TCanvas("c9","3D reconstruction of theta angle",75,75,300,1050); + c9->Divide(1,3); + + // data per hit + c9->cd(1); + hitsTheta->SetFillColor(5); + if (evNumber2>10) + hitsTheta->Fit("gaus"); + hitsTheta->Draw(); + + //data per track + c9->cd(2); + OriginalTheta->SetFillColor(5); + if (evNumber2>10) + OriginalTheta->Fit("gaus"); + OriginalTheta->Draw(); + + //recontructed data + c9->cd(3); Theta->SetFillColor(5); + if (evNumber2>10) + Theta->Fit("gaus"); Theta->Draw(); - c8->cd(7); + + c10 = new TCanvas("c10","3D reconstruction of cherenkov angle",100,100,300,1050); + c10->Divide(1,3); + + // data per hit + c10->cd(1); + ckovangle->SetFillColor(5); + ckovangle->SetXTitle("angle (radians)"); + if (evNumber2>10) + ckovangle->Fit("gaus"); + ckovangle->Draw(); + + //data per track + c10->cd(2); + OriginalOmega->SetFillColor(5); + OriginalOmega->SetXTitle("angle (radians)"); + if (evNumber2>10) + OriginalOmega->Fit("gaus"); + OriginalOmega->Draw(); + + //recontructed data + c10->cd(3); Omega3D->SetFillColor(5); Omega3D->SetXTitle("angle (radians)"); + if (evNumber2>10) + Omega3D->Fit("gaus"); Omega3D->Draw(); - c8->cd(8); + + + c11 = new TCanvas("c11","3D reconstruction of mean radius",125,125,300,700); + c11->Divide(1,2); + + // data per hit + c11->cd(1); + radius->SetFillColor(5); + radius->SetXTitle("radius (cm)"); + radius->Draw(); + + //recontructed data + c11->cd(2); MeanRadius->SetFillColor(5); MeanRadius->SetXTitle("radius (cm)"); MeanRadius->Draw(); + + + c12 = new TCanvas("c12","Cerenkov angle vs. Momentum",150,150,550,350); + + c12->cd(1); + identification->SetFillColor(5); + identification->SetXTitle("Momentum (GeV/c)"); + identification->SetYTitle("Cherenkov angle (radians)"); + + //Float_t pionmass=.139; + //Float_t kaonmass=.493; + //Float_t protonmass=.938; + //Float_t n=1.295; + + TF1 *pionplot = new TF1("pion","acos(sqrt((.139*.139+x*x)/(x*x*1.285*1.285)))",1,5); + TF1 *kaonplot = new TF1("kaon","acos(sqrt((.439*.439+x*x)/(x*x*1.285*1.285)))",1,5); + TF1 *protonplot = new TF1("proton","acos(sqrt((.938*.938+x*x)/(x*x*1.285*1.285)))",1,5); + + identification->Draw(); + + pionplot->SetLineColor(5); + pionplot->Draw("same"); + + kaonplot->SetLineColor(4); + kaonplot->Draw("same"); + + protonplot->SetLineColor(3); + protonplot->Draw("same"); + //identification->Draw("same"); + + + + c13 = new TCanvas("c13","Reconstruction Errors",200,200,900,350); + c13->Divide(3,1); + + c13->cd(1); + PhiError->SetFillColor(5); + if (evNumber2>10) + PhiError->Fit("gaus"); + PhiError->Draw(); + c13->cd(2); + ThetaError->SetFillColor(5); + if (evNumber2>10) + ThetaError->Fit("gaus"); + ThetaError->Draw(); + c13->cd(3); + OmegaError->SetFillColor(5); + OmegaError->SetXTitle("angle (radians)"); + if (evNumber2>10) + OmegaError->Fit("gaus"); + OmegaError->Draw(); } @@ -4075,3 +4419,27 @@ AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); printf("\nEnd of analysis\n"); printf("**********************************\n"); } + +//////////////////////////////////////////////////////////////////////// +void AliRICH::MakeBranchInTreeD(TTree *treeD, const char *file) +{ + // + // Create TreeD branches for the RICH. + // + + const Int_t kBufferSize = 4000; + char branchname[30]; + + // + // one branch for digits per chamber + // + for (Int_t i=0; i