X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RICH%2FAliRICH.cxx;h=924506f3031faeaa11460c12ad10567ab2b61cb6;hb=f701e8863eb28607018eb7ae3a8f151bee52798b;hp=05bb86d275202f0c3760c8bdfa4aeca2cd26739b;hpb=6197ac87943502b45732bce730ea62dc5baa3f73;p=u%2Fmrichter%2FAliRoot.git diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index 05bb86d2752..924506f3031 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -15,6 +15,89 @@ /* $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) + + Revision 1.47 2001/03/14 18:13:56 jbarbosa + Several changes to adapt to new IO. + Removed digitising function, using AliRICHMerger::Digitise from now on. + + Revision 1.46 2001/03/12 17:46:33 hristov + Changes needed on Sun with CC 5.0 + + Revision 1.45 2001/02/27 22:11:46 jbarbosa + Testing TreeS, removing of output. + + Revision 1.44 2001/02/27 15:19:12 jbarbosa + Transition to SDigits. + + Revision 1.43 2001/02/23 17:19:06 jbarbosa + Corrected photocathode definition in BuildGeometry(). + + Revision 1.42 2001/02/13 20:07:23 jbarbosa + Parametrised definition of photcathode dimensions. New spacers. New data members in AliRICHHit to store particle momentum + when entering the freon. Corrected calls to particle stack. + + Revision 1.41 2001/01/26 20:00:20 hristov + Major upgrade of AliRoot code + + Revision 1.40 2001/01/24 20:58:03 jbarbosa + Enhanced BuildGeometry. Now the photocathodes are drawn. + + Revision 1.39 2001/01/22 21:40:24 jbarbosa + Removing magic numbers + + Revision 1.37 2000/12/20 14:07:25 jbarbosa + Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova) + + Revision 1.36 2000/12/18 17:45:54 jbarbosa + Cleaned up PadHits object. + + Revision 1.35 2000/12/15 16:49:40 jbarbosa + Geometry and materials updates (wire supports, pcbs, backplane supports, frame). + Revision 1.34 2000/11/10 18:12:12 jbarbosa Bug fix for AliRICHCerenkov (thanks to P. Hristov) @@ -96,15 +179,21 @@ #include #include #include +#include +#include +#include +#include +#include -#include +#include #include #include "AliRICH.h" #include "AliSegmentation.h" +#include "AliRICHSegmentationV0.h" #include "AliRICHHit.h" #include "AliRICHCerenkov.h" -#include "AliRICHPadHit.h" +#include "AliRICHSDigit.h" #include "AliRICHDigit.h" #include "AliRICHTransientDigit.h" #include "AliRICHRawCluster.h" @@ -112,34 +201,41 @@ #include "AliRICHRecHit3D.h" #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" -#include "TGeant3.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; -static TClonesArray *fClusters2; -static TClonesArray *fHits2; -static TTree *TrH1; 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; - fPadHits = 0; - fNPadHits = 0; + fSDigits = 0; + fNSDigits = 0; fNcerenkovs = 0; fDchambers = 0; fRecHits1D = 0; @@ -147,34 +243,29 @@ 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); - fPadHits = new TClonesArray("AliRICHPadHit",100000); + fSDigits = new TClonesArray("AliRICHSDigit",100000); fCerenkovs = new TClonesArray("AliRICHCerenkov",1000); gAlice->AddHitList(fCerenkovs); - //gAlice->AddHitList(fHits); - fNPadHits = 0; + fNSDigits = 0; fNcerenkovs = 0; fIshunt = 0; @@ -188,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; } @@ -198,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); @@ -220,81 +315,209 @@ 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; - delete fPadHits; + delete fSDigits; delete fCerenkovs; + + //PH Delete TObjArrays + if (fChambers) { + fChambers->Delete(); + delete fChambers; + } + if (fDchambers) { + fDchambers->Delete(); + delete fDchambers; + } + if (fRawClusters) { + fRawClusters->Delete(); + delete fRawClusters; + } + if (fRecHits1D) { + fRecHits1D->Delete(); + delete fRecHits1D; + } + if (fRecHits3D) { + fRecHits3D->Delete(); + delete fRecHits3D; + } + } -//___________________________________________ -void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits) + +//_____________________________________________________________________________ +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 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; -// -// Adds a hit to the Hits list + ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, nnew, newclust, res); + Int_t ic=0; + // +// Add new clusters + for (Int_t i=0; i 0) { + ic++; +// Cluster Charge + clhits[1] = Int_t(newclust[0][i]); +// Pad: ix + clhits[2] = Int_t(newclust[1][i]); +// Pad: iy + clhits[3] = Int_t(newclust[2][i]); +// Pad: chamber sector + clhits[4] = Int_t(newclust[3][i]); + + //printf(" %d %d %d %d %d\n", clhits[0], clhits[1], clhits[2], clhits[3], clhits[4]); + + AddSDigit(clhits); + } + } + + if (gAlice->TreeS()){ + gAlice->TreeS()->Fill(); + gAlice->TreeS()->Write(0,TObject::kOverwrite); + //printf("Filled SDigits...\n"); + } + + 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. + if(IsDebugHit()||IsDebugDigit()) cout<\n"; + + + int nparticles = gAlice->GetNtrack(); + cout << "Particles (RICH):" < 0) printf("SDigits were already generated.\n"); + +} + +//___________________________________________ +void AliRICH::SDigits2Digits(Int_t nev, Int_t flag) +{ +// Generate digits. +// Called from macro. Multiple events, more functionality. + if(IsDebugDigit()) cout<\n"; + + //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):" <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() +{ + SDigits2Digits(0,0); +} +//___________________________________________ +void AliRICH::Digits2Reco() +{ +// Generate clusters +// Called from alirun, single event only. + if(IsDebugDigit()||IsDebugReco()) cout<\n"; + + + int nparticles = gAlice->GetNtrack(); + cout << "Particles (RICH):" < 0) FindClusters(0,0); + +} + +void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits) +{ +// Adds the current hit to the RICH hits list + if(IsDebugHit()) cout<\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::AddPadHit(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"; - TClonesArray &lPadHits = *fPadHits; - new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits); + TClonesArray &lSDigits = *fSDigits; + 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 - // - - 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); } @@ -306,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); } //___________________________________________ @@ -330,75 +551,251 @@ void AliRICH::BuildGeometry() // // Builds a TNode geometry for event display // - TNode *node, *top; + TNode *node, *subnode, *top; - const int kColorRICH = kGreen; + const int kColorRICH = kRed; // top=gAlice->GetGeometry()->GetNode("alice"); - + + AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH"); + AliRICHSegmentationV0* segmentation; + AliRICHChamber* iChamber; + AliRICHGeometry* geometry; + + iChamber = &(pRICH->Chamber(0)); + segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(); + geometry=iChamber->GetGeometryModel(); new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15); - + + Float_t padplane_width = segmentation->GetPadPlaneWidth(); + Float_t padplane_length = segmentation->GetPadPlaneLength(); + + //printf("\n\n\n\n\n In BuildGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f\n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy()); + + new TBRIK("PHOTO","PHOTO","void", padplane_width/2,.1,padplane_length/2); + + //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", padplane_width/2,padplane_length/2); + //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", segmentation->GetPadPlaneWidth(), segmentation->GetPadPlaneLength()); + + Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane + Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction + Float_t deltatheta = 20; //theta angle between center of chambers - x direction + Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180); + Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180); + Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180); + Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180); + + //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta); + + new TRotMatrix("rot993","rot993",90., 0. , 90. - deltaphi, 90. , deltaphi, -90. ); + new TRotMatrix("rot994","rot994",90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. ); + new TRotMatrix("rot995","rot995",90., 0. , 90. , 90. , 0. , 0. ); + new TRotMatrix("rot996","rot996",90., deltatheta , 90. , 90 + deltatheta , 0. , 0. ); + new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta); + new TRotMatrix("rot998","rot998",90., 0. , 90 + deltaphi , 90. , deltaphi, 90. ); + new TRotMatrix("rot999","rot999",90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta); + + Float_t pos1[3]={0. , offset*cosphi , offset*sinphi}; + Float_t pos2[3]={offset*sintheta , offset*costheta , 0. }; + Float_t pos3[3]={0. , offset , 0.}; + Float_t pos4[3]={-offset*sintheta , offset*costheta , 0.}; + Float_t pos5[3]={offset*sinphi , offset*costheta*cosphi, -offset*sinphi}; + Float_t pos6[3]={0. , offset*cosphi , -offset*sinphi}; + Float_t pos7[3]={ -offset*sinphi , offset*costheta*cosphi, -offset*sinphi}; + + top->cd(); - Float_t pos1[3]={0,471.8999,165.2599}; + //Float_t pos1[3]={0,471.8999,165.2599}; //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2], - new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90); + //new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90); node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993"); - - node->SetLineColor(kColorRICH); + node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); fNodes->Add(node); - top->cd(); - - Float_t pos2[3]={171,470,0}; + + + top->cd(); + //Float_t pos2[3]={171,470,0}; //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2], - new TRotMatrix("rot994","rot994",90,-20,90,70,0,0); + //new TRotMatrix("rot994","rot994",90,-20,90,70,0,0); node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994"); - - node->SetLineColor(kColorRICH); + node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); fNodes->Add(node); + + top->cd(); - Float_t pos3[3]={0,500,0}; + //Float_t pos3[3]={0,500,0}; //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2], - new TRotMatrix("rot995","rot995",90,0,90,90,0,0); + //new TRotMatrix("rot995","rot995",90,0,90,90,0,0); node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995"); - - node->SetLineColor(kColorRICH); + node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); fNodes->Add(node); + top->cd(); - Float_t pos4[3]={-171,470,0}; + //Float_t pos4[3]={-171,470,0}; //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2], - new TRotMatrix("rot996","rot996",90,20,90,110,0,0); + //new TRotMatrix("rot996","rot996",90,20,90,110,0,0); node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996"); - - node->SetLineColor(kColorRICH); + node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); fNodes->Add(node); + + top->cd(); - Float_t pos5[3]={161.3999,443.3999,-165.3}; + //Float_t pos5[3]={161.3999,443.3999,-165.3}; //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2], - new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70); + //new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70); node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997"); - node->SetLineColor(kColorRICH); + node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); fNodes->Add(node); + + top->cd(); - Float_t pos6[3]={0., 471.9, -165.3,}; + //Float_t pos6[3]={0., 471.9, -165.3,}; //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2], - new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90); + //new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90); node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998"); - - node->SetLineColor(kColorRICH); - fNodes->Add(node); + fNodes->Add(node);node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + + top->cd(); - Float_t pos7[3]={-161.399,443.3999,-165.3}; + //Float_t pos7[3]={-161.399,443.3999,-165.3}; //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2], - new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110); + //new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110); node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999"); node->SetLineColor(kColorRICH); + node->cd(); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); + subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,""); + subnode->SetLineColor(kGreen); + fNodes->Add(subnode); fNodes->Add(node); } @@ -425,12 +822,12 @@ void AliRICH::CreateGeometry() //End_Html AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH"); - AliSegmentation* segmentation; + AliRICHSegmentationV0* segmentation; AliRICHGeometry* geometry; AliRICHChamber* iChamber; iChamber = &(pRICH->Chamber(0)); - segmentation=iChamber->GetSegmentationModel(0); + segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(); geometry=iChamber->GetGeometryModel(); Float_t distance; @@ -441,10 +838,15 @@ void AliRICH::CreateGeometry() Float_t oqua_thickness = .5; //CsI dimensions - Float_t csi_length = 160*.8 + 2.6; - Float_t csi_width = 144*.84 + 2*2.6; - - Int_t *idtmed = fIdtmed->GetArray()-999; + //Float_t csi_length = 160*.8 + 2.6; + //Float_t csi_width = 144*.84 + 2*2.6; + + Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone(); + Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone(); + + //printf("\n\n\n\n\n In CreateGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f deadzone: %f \n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy(),segmentation->DeadZone()); + + Int_t *idtmed = fIdtmed->GetArray()-999; Int_t i; Float_t zs; @@ -729,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"); @@ -767,58 +1169,48 @@ void AliRICH::CreateGeometry() AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.); - //Int_t nspacers = (Int_t)(TMath::Abs(geometry->GetInnerFreonLength()/14.4)); - Int_t nspacers = 9; + //Placing of the spacers inside the freon slabs + + Int_t nspacers = 30; //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers); //printf("Nspacers: %d", nspacers); - //for (i = 1; i <= 9; ++i) { - //zs = (5 - i) * 14.4; //Original settings - for (i = 0; i < nspacers; i++) { - zs = (TMath::Abs(nspacers/2) - i) * 14.4; - gMC->Gspos("SPAC", i, "FRE1", 6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings - //gMC->Gspos("SPAC", i, "FRE1", zs, 0., 6.7, idrotm[1019], "ONLY"); + for (i = 0; i < nspacers/3; i++) { + zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2; + gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings + } + + for (i = nspacers/3; i < (nspacers*2)/3; i++) { + zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2; + gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings } - //for (i = 10; i <= 18; ++i) { - //zs = (14 - i) * 14.4; //Original settings - for (i = nspacers; i < nspacers*2; ++i) { - zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4; - gMC->Gspos("SPAC", i, "FRE1", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings - //gMC->Gspos("SPAC", i, "FRE1", zs, 0., -6.7, idrotm[1019], "ONLY"); + + for (i = (nspacers*2)/3; i < nspacers; ++i) { + zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2; + gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings } - //for (i = 1; i <= 9; ++i) { - //zs = (5 - i) * 14.4; //Original settings - for (i = 0; i < nspacers; i++) { - zs = (TMath::Abs(nspacers/2) - i) * 14.4; - gMC->Gspos("SPAC", i, "FRE2", 6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings - //gMC->Gspos("SPAC", i, "FRE2", zs, 0., 6.7, idrotm[1019], "ONLY"); + for (i = 0; i < nspacers/3; i++) { + zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2; + gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings } - //for (i = 10; i <= 18; ++i) { - //zs = (5 - i) * 14.4; //Original settings - for (i = nspacers; i < nspacers*2; ++i) { - zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4; - gMC->Gspos("SPAC", i, "FRE2", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings - //gMC->Gspos("SPAC", i, "FRE2", zs, 0., -6.7, idrotm[1019], "ONLY"); + + for (i = nspacers/3; i < (nspacers*2)/3; i++) { + zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2; + gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings } - /*gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY"); - gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY"); - gMC->Gspos("OQF1", 1, "SRIC", 31.3, -4.724, 41.3, 0, "ONLY"); - gMC->Gspos("OQF2", 2, "SRIC", 0., -4.724, 0., 0, "ONLY"); - gMC->Gspos("OQF1", 3, "SRIC", -31.3, -4.724, -41.3, 0, "ONLY"); - gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY"); //Original settings - gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY"); //Original settings - gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY"); - gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY"); - gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY"); - gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");*/ - + for (i = (nspacers*2)/3; i < nspacers; ++i) { + zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2; + gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings + } + gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY"); gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY"); gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2 + 2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3) +// printf("Opaque quartz in SRIC %f\n", 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2); gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2) - 2, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (-31.3) //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (-21.65) @@ -827,6 +1219,7 @@ void AliRICH::CreateGeometry() gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY"); gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY"); gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY"); + printf("CSI pos: %f\n",1.276 + geometry->GetGapThickness()/2 + .25); // Wire support placing @@ -845,30 +1238,60 @@ 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"); //printf("Position of the gap: %f to %f\n", 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - .2, 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 + .2); // Place RICH inside ALICE apparatus - - AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.); - AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.); - AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.); - AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.); - AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.); - AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.); - AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.); - - gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY"); - gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY"); - gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY"); - gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY"); - gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY"); - gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY"); - gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY"); + + /* old values + + AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.); + AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.); + AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.); + AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.); + AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.); + AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.); + AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.); + + gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY"); + gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY"); + gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY"); + gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY"); + gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY"); + gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY"); + gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");*/ + + // The placing of the chambers is measured from the vertex to the base of the methane vessel (490 cm) + + Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane + Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction + Float_t deltatheta = 20; //theta angle between center of chambers - x direction + Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180); + Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180); + Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180); + Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180); + + //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta); + + AliMatrix(idrotm[1000], 90., 0. , 90. - deltaphi, 90. , deltaphi, -90. ); + AliMatrix(idrotm[1001], 90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. ); + AliMatrix(idrotm[1002], 90., 0. , 90. , 90. , 0. , 0. ); + AliMatrix(idrotm[1003], 90., deltatheta , 90. , 90 + deltatheta , 0. , 0. ); + AliMatrix(idrotm[1004], 90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta); + AliMatrix(idrotm[1005], 90., 0. , 90 + deltaphi , 90. , deltaphi, 90. ); + AliMatrix(idrotm[1006], 90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta); + + gMC->Gspos("RICH", 1, "ALIC", 0. , offset*cosphi , offset*sinphi ,idrotm[1000], "ONLY"); + gMC->Gspos("RICH", 2, "ALIC", (offset)*sintheta , offset*costheta , 0. ,idrotm[1001], "ONLY"); + gMC->Gspos("RICH", 3, "ALIC", 0. , offset , 0. ,idrotm[1002], "ONLY"); + gMC->Gspos("RICH", 4, "ALIC", -(offset)*sintheta, offset*costheta , 0. ,idrotm[1003], "ONLY"); + gMC->Gspos("RICH", 5, "ALIC", (offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1004], "ONLY"); + gMC->Gspos("RICH", 6, "ALIC", 0. , offset*cosphi , -offset*sinphi,idrotm[1005], "ONLY"); + gMC->Gspos("RICH", 7, "ALIC", -(offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1006], "ONLY"); } @@ -1068,8 +1491,6 @@ void AliRICH::CreateMaterials() Int_t *idtmed = fIdtmed->GetArray()-999; - TGeant3 *geant3 = (TGeant3*) gMC; - // --- Photon energy (GeV) // --- Refraction indexes for (i = 0; i < 26; ++i) { @@ -1186,17 +1607,17 @@ void AliRICH::CreateMaterials() AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - geant3->Gsckov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane); - geant3->Gsckov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane); - geant3->Gsckov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz); - geant3->Gsckov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon); - geant3->Gsckov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane); - geant3->Gsckov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane); - geant3->Gsckov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid); - geant3->Gsckov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz); - geant3->Gsckov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane); - geant3->Gsckov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid); - geant3->Gsckov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz); + gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane); + gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane); + gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz); + gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon); + gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane); + gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane); + gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid); + gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz); + gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane); + gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid); + gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz); } //___________________________________________ @@ -1417,67 +1838,100 @@ Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t ) } //___________________________________________ -void AliRICH::MakeBranch(Option_t* option) +void AliRICH::MakeBranch(Option_t* option, const char *file) { // Create Tree branches for the RICH. const Int_t kBufferSize = 4000; char branchname[20]; - - - AliDetector::MakeBranch(option); - sprintf(branchname,"%sCerenkov",GetName()); - if (fCerenkovs && gAlice->TreeH()) { - gAlice->TreeH()->Branch(branchname,&fCerenkovs, kBufferSize); - printf("Making Branch %s for Cerenkov Hits\n",branchname); - } - - sprintf(branchname,"%sPadHits",GetName()); - if (fPadHits && gAlice->TreeH()) { - gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize); - printf("Making Branch %s for PadHits\n",branchname); + + AliDetector::MakeBranch(option,file); + + const char *cH = strstr(option,"H"); + const char *cD = strstr(option,"D"); + const char *cR = strstr(option,"R"); + const char *cS = strstr(option,"S"); + + + if (cH) { + sprintf(branchname,"%sCerenkov",GetName()); + if (fCerenkovs && gAlice->TreeH()) { + //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 = 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()); + } + } + + if (cS) { + sprintf(branchname,"%sSDigits",GetName()); + if (fSDigits && gAlice->TreeS()) { + //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()); + } } -// one branch for digits per chamber - Int_t i; + if (cD) { + // + // one branch for digits per chamber + // + Int_t i; - for (i=0; iTreeD()) { - gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize); - printf("Making Branch %s for digits in chamber %d\n",branchname,i+1); + //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); } + } } -// one branch for raw clusters per chamber - for (i=0; iTreeR()) { - gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize); - printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1); - } - } + if (cR) { + // + // one branch for raw clusters per chamber + // - // one branch for rec hits per chamber - for (i=0; iTreeR()) { - gAlice->TreeR()->Branch(branchname,&((*fRecHits1D)[i]), kBufferSize); - printf("Making Branch %s for 1D rec. hits in chamber %d\n",branchname,i+1); - } - } - for (i=0; iTreeR()) { - gAlice->TreeR()->Branch(branchname,&((*fRecHits3D)[i]), kBufferSize); - printf("Making Branch %s for 3D rec. hits in chamber %d\n",branchname,i+1); - } - } - + //printf("Called MakeBranch for TreeR\n"); + + Int_t i; + + for (i=0; iTreeR()) { + //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ; + MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ; + //branch->SetAutoDelete(kFALSE); + } + } + // + // one branch for rec hits per chamber + // + for (i=0; iTreeR()) { + //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()) { + MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ; + //branch->SetAutoDelete(kFALSE); + } + } + } } //___________________________________________ @@ -1493,24 +1947,41 @@ void AliRICH::SetTreeAddress() TTree *treeH = gAlice->TreeH(); TTree *treeD = gAlice->TreeD(); TTree *treeR = gAlice->TreeR(); + TTree *treeS = gAlice->TreeS(); if (treeH) { - if (fPadHits) { - branch = treeH->GetBranch("RICHPadHits"); - if (branch) branch->SetAddress(&fPadHits); - } - if (fCerenkovs) { + if (fCerenkovs) { branch = treeH->GetBranch("RICHCerenkov"); if (branch) branch->SetAddress(&fCerenkovs); } + if (fSDigits) { + branch = treeH->GetBranch("RICHSDigits"); + if (branch) + { + branch->SetAddress(&fSDigits); + //printf("Setting sdigits branch address at %p in TreeH\n",&fSDigits); + } + } + } + + if (treeS) { + if (fSDigits) { + branch = treeS->GetBranch("RICHSDigits"); + if (branch) + { + branch->SetAddress(&fSDigits); + //printf("Setting sdigits branch address at %p in TreeS\n",&fSDigits); + } + } } + if (treeD) { for (int i=0; iGetBranch(branchname); - if (branch) branch->SetAddress(&((*fDchambers)[i])); + branch = treeD->GetBranch(branchname); + if (branch) branch->SetAddress(&((*fDchambers)[i])); } } } @@ -1546,9 +2017,9 @@ void AliRICH::ResetHits() { // Reset number of clusters and the cluster array for this detector AliDetector::ResetHits(); - fNPadHits = 0; + fNSDigits = 0; fNcerenkovs = 0; - if (fPadHits) fPadHits->Clear(); + if (fSDigits) fSDigits->Clear(); if (fCerenkovs) fCerenkovs->Clear(); } @@ -1560,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; } } @@ -1572,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; } } @@ -1585,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; } } @@ -1598,77 +2072,23 @@ 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) +void AliRICH::StepManager() { - -// -// 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::SetNsec(Int_t id, Int_t nsec) -{ - -// -// Sets the number of padplanes -// - - ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec); -} - - -//___________________________________________ -void AliRICH::StepManager() -{ - -// Full Step Manager +// Full Step Manager Int_t copy, id; static Int_t idvol; static Int_t vol[2]; Int_t ipart; - static Float_t hits[18]; + static Float_t hits[22]; static Float_t ckovData[19]; TLorentzVector position; TLorentzVector momentum; @@ -1679,14 +2099,13 @@ 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; const Float_t kBig=1.e10; TClonesArray &lhits = *fHits; - TGeant3 *geant3 = (TGeant3*) gMC; TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()]; //if (current->Energy()>1) @@ -1730,7 +2149,7 @@ void AliRICH::StepManager() //printf("Track entered (1)\n"); if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) { //is it in freo? - if (geant3->Gctrak()->nstep<1){ //is it the first step? + if (gMC->IsNewTrack()){ //is it the first step? //printf("I'm in!\n"); Int_t mother = current->GetFirstMother(); @@ -1746,7 +2165,7 @@ void AliRICH::StepManager() } //first step question } //freo question - if (geant3->Gctrak()->nstep<1){ //is it first step? + if (gMC->IsNewTrack()){ //is it first step? if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz? { ckovData[12] = 2; @@ -1779,10 +2198,11 @@ 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) { - geant3->StopTrack(); + gMC->StopTrack(); ckovData[13] = 5; AddCerenkov(gAlice->CurrentTrack(),vol,ckovData); //printf("Added One (1)!\n"); @@ -1807,9 +2227,10 @@ 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) { - geant3->StopTrack(); + gMC->StopTrack(); ckovData[13] = 6; AddCerenkov(gAlice->CurrentTrack(),vol,ckovData); //printf("Added One (2)!\n"); @@ -1823,21 +2244,22 @@ void AliRICH::StepManager() /********************Evaluation of losses************************/ /******************still in the old fashion**********************/ - Int_t i1 = geant3->Gctrak()->nmec; //number of physics mechanisms acting on the particle + TArrayI procs; + Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle for (Int_t i = 0; i < i1; ++i) { // Reflection loss - if (geant3->Gctrak()->lmec[i] == 106) { //was it reflected + if (procs[i] == kPLightReflection) { //was it reflected ckovData[13]=10; if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) ckovData[13]=1; if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) ckovData[13]=2; - //geant3->StopTrack(); + //gMC->StopTrack(); //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData); } //reflection question // Absorption loss - else if (geant3->Gctrak()->lmec[i] == 101) { //was it absorbed? + else if (procs[i] == kPLightAbsorption) { //was it absorbed? //printf("Got in absorption\n"); ckovData[13]=20; if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) @@ -1856,7 +2278,7 @@ void AliRICH::StepManager() if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) { ckovData[13]=16; } - geant3->StopTrack(); + gMC->StopTrack(); AddCerenkov(gAlice->CurrentTrack(),vol,ckovData); //printf("Added One (3)!\n"); //printf("Added cerenkov %d\n",fCkovNumber); @@ -1864,9 +2286,9 @@ void AliRICH::StepManager() // Photon goes out of tracking scope - else if (geant3->Gctrak()->lmec[i] == 30) { //is it below energy treshold? + else if (procs[i] == kPStop) { //is it below energy treshold? ckovData[13]=21; - geant3->StopTrack(); + gMC->StopTrack(); AddCerenkov(gAlice->CurrentTrack(),vol,ckovData); //printf("Added One (4)!\n"); } // energy treshold question @@ -1884,11 +2306,16 @@ void AliRICH::StepManager() if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) { //printf("Cerenkov\n"); - if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) + + //if (gMC->TrackPid() == 50000051) + //printf("Tracking a feedback\n"); + + if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) { //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName()); //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy)); //printf("Got in CSI\n"); + //printf("Tracking a %d\n",gMC->TrackPid()); if (gMC->Edep() > 0.){ gMC->TrackPosition(position); gMC->TrackMomentum(momentum); @@ -1903,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; @@ -1919,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 @@ -1928,7 +2367,7 @@ void AliRICH::StepManager() ckovData[3] = pos[2]; // Z-position for hit ckovData[4] = theta; // theta angle of incidence ckovData[5] = phi; // phi angle of incidence - ckovData[8] = (Float_t) fNPadHits; // first padhit + ckovData[8] = (Float_t) fNSDigits; // first sdigit ckovData[9] = -1; // last pad hit ckovData[13] = 4; // photon was detected ckovData[14] = mom[0]; @@ -1940,12 +2379,15 @@ void AliRICH::StepManager() cherenkovLoss += destep; ckovData[7]=cherenkovLoss; - nPads = MakePadHits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov); - if (fNPadHits > (Int_t)ckovData[8]) { + nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov); + + if (fNSDigits > (Int_t)ckovData[8]) { ckovData[8]= ckovData[8]+1; - ckovData[9]= (Float_t) fNPadHits; + ckovData[9]= (Float_t) fNSDigits; } + //printf("Cerenkov loss: %f\n", cherenkovLoss); + ckovData[17] = nPads; //printf("nPads:%d",nPads); @@ -1956,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); @@ -2001,6 +2443,14 @@ void AliRICH::StepManager() }*/ if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) { + gMC->TrackMomentum(momentum); + mom[0]=momentum(0); + mom[1]=momentum(1); + mom[2]=momentum(2); + mom[3]=momentum(3); + hits [19] = mom[0]; + hits [20] = mom[1]; + hits [21] = mom[2]; fFreonProd=1; } @@ -2024,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(); // @@ -2054,12 +2510,13 @@ void AliRICH::StepManager() hits[3] = localPos[2]; // Z-position for hit hits[4] = localTheta; // theta angle of incidence hits[5] = localPhi; // phi angle of incidence - hits[8] = (Float_t) fNPadHits; // first padhit + hits[8] = (Float_t) fNSDigits; // first sdigit hits[9] = -1; // last pad hit hits[13] = fFreonProd; // did id hit the freon? hits[14] = mom[0]; hits[15] = mom[1]; hits[16] = mom[2]; + hits[18] = 0; // dummy cerenkov angle tlength = 0; eloss = 0; @@ -2075,7 +2532,8 @@ void AliRICH::StepManager() if(idvolAt(idvol)) ->SigGenInit(localPos[0], localPos[2], localPos[1]); } } @@ -2096,7 +2554,7 @@ void AliRICH::StepManager() { if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n"); - nPads = MakePadHits(xhit,yhit,eloss,idvol,kMip); + nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); hits[17] = nPads; //printf("nPads:%d",nPads); } @@ -2104,9 +2562,9 @@ void AliRICH::StepManager() hits[6]=tlength; hits[7]=eloss; - if (fNPadHits > (Int_t)hits[8]) { + if (fNSDigits > (Int_t)hits[8]) { hits[8]= hits[8]+1; - hits[9]= (Float_t) fNPadHits; + hits[9]= (Float_t) fNSDigits; } //if(sector !=-1) @@ -2117,16 +2575,18 @@ 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) { if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n"); - nPads = MakePadHits(xhit,yhit,eloss,idvol,kMip); + nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); hits[17] = nPads; //printf("Npads:%d",NPads); } @@ -2144,7 +2604,7 @@ void AliRICH::StepManager() } /*************************************************End of MIP treatment**************************************/ //} -} +}//void AliRICH::StepManager() void AliRICH::FindClusters(Int_t nev,Int_t lastEntry) { @@ -2154,9 +2614,10 @@ void AliRICH::FindClusters(Int_t nev,Int_t lastEntry) // for (Int_t icat=1;icat<2;icat++) { gAlice->ResetDigits(); - gAlice->TreeD()->GetEvent(1); // spurious +1 ... + gAlice->TreeD()->GetEvent(0); for (Int_t ich=0;ichAt(ich); TClonesArray *pRICHdigits = this->DigitsAddress(ich); if (pRICHdigits == 0) continue; @@ -2199,625 +2660,1786 @@ void AliRICH::FindClusters(Int_t nev,Int_t lastEntry) //gObjectTable->Print(); } - -//______________________________________________________________________________ -void AliRICH::Streamer(TBuffer &R__b) -{ - // Stream an object of class AliRICH. - AliRICHChamber *iChamber; - AliSegmentation *segmentation; - AliRICHResponse *response; - TClonesArray *digitsaddress; - TClonesArray *rawcladdress; - TClonesArray *rechitaddress1D; - TClonesArray *rechitaddress3D; - - if (R__b.IsReading()) { - Version_t R__v = R__b.ReadVersion(); if (R__v) { } - AliDetector::Streamer(R__b); - R__b >> fNPadHits; - R__b >> fPadHits; // diff - R__b >> fNcerenkovs; - R__b >> fCerenkovs; // diff - R__b >> fDchambers; - R__b >> fRawClusters; - R__b >> fRecHits1D; //diff - R__b >> fRecHits3D; //diff - R__b >> fDebugLevel; //diff - R__b.ReadStaticArray(fNdch); - R__b.ReadStaticArray(fNrawch); - R__b.ReadStaticArray(fNrechits1D); - R__b.ReadStaticArray(fNrechits3D); -// - R__b >> fChambers; -// Stream chamber related information - for (Int_t i =0; iStreamer(R__b); - segmentation=iChamber->GetSegmentationModel(); - segmentation->Streamer(R__b); - response=iChamber->GetResponseModel(); - response->Streamer(R__b); - rawcladdress=(TClonesArray*) (*fRawClusters)[i]; - rawcladdress->Streamer(R__b); - rechitaddress1D=(TClonesArray*) (*fRecHits1D)[i]; - rechitaddress1D->Streamer(R__b); - rechitaddress3D=(TClonesArray*) (*fRecHits3D)[i]; - rechitaddress3D->Streamer(R__b); - digitsaddress=(TClonesArray*) (*fDchambers)[i]; - digitsaddress->Streamer(R__b); - } - R__b >> fDebugLevel; - R__b >> fCkovNumber; - R__b >> fCkovQuarz; - R__b >> fCkovGap; - R__b >> fCkovCsi; - R__b >> fLostRfreo; - R__b >> fLostRquar; - R__b >> fLostAfreo; - R__b >> fLostAquarz; - R__b >> fLostAmeta; - R__b >> fLostCsi; - R__b >> fLostWires; - R__b >> fFreonProd; - R__b >> fMipx; - R__b >> fMipy; - R__b >> fFeedbacks; - R__b >> fLostFresnel; - - } else { - R__b.WriteVersion(AliRICH::IsA()); - AliDetector::Streamer(R__b); - R__b << fNPadHits; - R__b << fPadHits; // diff - R__b << fNcerenkovs; - R__b << fCerenkovs; // diff - R__b << fDchambers; - R__b << fRawClusters; - R__b << fRecHits1D; //diff - R__b << fRecHits3D; //diff - R__b << fDebugLevel; //diff - R__b.WriteArray(fNdch, kNCH); - R__b.WriteArray(fNrawch, kNCH); - R__b.WriteArray(fNrechits1D, kNCH); - R__b.WriteArray(fNrechits3D, kNCH); -// - R__b << fChambers; -// Stream chamber related information - for (Int_t i =0; iStreamer(R__b); - segmentation=iChamber->GetSegmentationModel(); - segmentation->Streamer(R__b); - response=iChamber->GetResponseModel(); - response->Streamer(R__b); - rawcladdress=(TClonesArray*) (*fRawClusters)[i]; - rawcladdress->Streamer(R__b); - rechitaddress1D=(TClonesArray*) (*fRecHits1D)[i]; - rechitaddress1D->Streamer(R__b); - rechitaddress3D=(TClonesArray*) (*fRecHits3D)[i]; - rechitaddress3D->Streamer(R__b); - digitsaddress=(TClonesArray*) (*fDchambers)[i]; - digitsaddress->Streamer(R__b); - } - R__b << fDebugLevel; - R__b << fCkovNumber; - R__b << fCkovQuarz; - R__b << fCkovGap; - R__b << fCkovCsi; - R__b << fLostRfreo; - R__b << fLostRquar; - R__b << fLostAfreo; - R__b << fLostAquarz; - R__b << fLostAmeta; - R__b << fLostCsi; - R__b << fLostWires; - R__b << fFreonProd; - R__b << fMipx; - R__b << fMipy; - R__b << fFeedbacks; - R__b << fLostFresnel; - } -} -AliRICHPadHit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters ) +AliRICHSDigit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters ) { // // Initialise the pad iterator - // Return the address of the first padhit for hit + // 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); - return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1); + if (nclust && hit->PHlast() > 0) { + sMaxIterPad=Int_t(hit->PHlast()); + sCurIterPad=Int_t(hit->PHfirst()); + return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1); } else { return 0; } } -AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters) +AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters) { // Iterates over pads sCurIterPad++; if (sCurIterPad <= sMaxIterPad) { - return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1); + return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1); } else { return 0; } } +AliRICH& AliRICH::operator=(const AliRICH& rhs) +{ +// Assignment operator + return *this; + +} -void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename) +void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2) { - // keep galice.root for signal and name differently the file for - // background when add! otherwise the track info for signal will be lost ! + + Int_t NpadX = 162; // number of pads on X + Int_t NpadY = 162; // number of pads on Y + + Int_t Pad[162][162]; + for (Int_t i=0;i3 GeV primary tracks",100,0,50); + TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600); + + + - static Bool_t first=kTRUE; - static TFile *pFile; - char *addBackground = strstr(option,"Add"); - Int_t particle; +// Start loop over events - FILE* points; //these will be the digits... + Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0; + Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0; + Int_t photons=0, primaryphotons=0, highprimaryphotons=0; + TRandom* random=0; - points=fopen("points.dat","w"); + for (int nev=0; nev<= evNumber2; nev++) { + Int_t nparticles = gAlice->GetEvent(nev); + - AliRICHChamber* iChamber; - AliSegmentation* segmentation; - - Int_t digitise=0; - Int_t trk[50]; - Int_t chtrk[50]; - TObjArray *list=new TObjArray; - static TClonesArray *pAddress=0; - if(!pAddress) pAddress=new TClonesArray("TVector",1000); - Int_t digits[5]; - - AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH"); - AliHitMap* pHitMap[10]; - Int_t i; - for (i=0; i<10; i++) {pHitMap[i]=0;} - if (addBackground ) { - if(first) { - fFileName=filename; - cout<<"filename"<cd(); - // Get Hits Tree header from file - if(fHits2) fHits2->Clear(); - if(fClusters2) fClusters2->Clear(); - if(TrH1) delete TrH1; - TrH1=0; - - char treeName[20]; - sprintf(treeName,"TreeH%d",nev); - TrH1 = (TTree*)gDirectory->Get(treeName); - if (!TrH1) { - printf("ERROR: cannot find Hits Tree for event:%d\n",nev); - } - // Set branch addresses - TBranch *branch; - char branchname[20]; - sprintf(branchname,"%s",GetName()); - if (TrH1 && fHits2) { - branch = TrH1->GetBranch(branchname); - if (branch) branch->SetAddress(&fHits2); - } - if (TrH1 && fClusters2) { - branch = TrH1->GetBranch("RICHCluster"); - if (branch) branch->SetAddress(&fClusters2); - } - } - - AliHitMap* hm; - Int_t countadr=0; - Int_t counter=0; - for (i =0; iGetSegmentationModel(1); - pHitMap[i] = new AliRICHHitMapA1(segmentation, list); - } - // - // Loop over tracks - // - - TTree *treeH = gAlice->TreeH(); - Int_t ntracks =(Int_t) treeH->GetEntries(); - for (Int_t track=0; trackResetHits(); - treeH->GetEvent(track); - // - // Loop over hits - for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1); - mHit; - mHit=(AliRICHHit*)pRICH->NextHit()) - { - - Int_t nch = mHit->fChamber-1; // chamber number - Int_t index = mHit->Track(); - if (nch >kNCH) continue; - iChamber = &(pRICH->Chamber(nch)); - - TParticle *current = (TParticle*)(*gAlice->Particles())[index]; - - if (current->GetPdgCode() >= 50000050) - { - TParticle *motherofcurrent = (TParticle*)(*gAlice->Particles())[current->GetFirstMother()]; - particle = motherofcurrent->GetPdgCode(); - } - else - { - particle = current->GetPdgCode(); - } + printf ("Event number : %d\n",nev); + printf ("Number of particles: %d\n",nparticles); + if (nev < evNumber1) continue; + if (nparticles <= 0) return; + +// Get pointers to RICH detector and Hits containers + + AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH"); + + TTree *treeH = gAlice->TreeH(); + Int_t ntracks =(Int_t) treeH->GetEntries(); + +// Start loop on tracks in the hits containers + + for (Int_t track=0; trackResetHits(); + treeH->GetEvent(track); + + for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1); + mHit; + mHit=(AliRICHHit*)pRICH->NextHit()) + { + //Int_t nch = mHit->fChamber; // chamber number + //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->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->Particle()); + Float_t R; + Float_t PTfinal; + Float_t PTvertex; + + TParticle *current = gAlice->Particle(index); + + //Float_t energy=current->Energy(); - - //printf("Flag:%d\n",flag); - //printf("Track:%d\n",track); - //printf("Particle:%d\n",particle); - - digitise=1; - - if (flag == 1) - if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111) - digitise=0; - - if (flag == 2) - if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 - || TMath::Abs(particle)==311) - digitise=0; - - if (flag == 3 && TMath::Abs(particle)==2212) - digitise=0; - - if (flag == 4 && TMath::Abs(particle)==13) - digitise=0; - - if (flag == 5 && TMath::Abs(particle)==11) - digitise=0; - - if (flag == 6 && TMath::Abs(particle)==2112) - digitise=0; - - - //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise); - - - if (digitise) - { + R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy()); + PTfinal=TMath::Sqrt(px*px + py*py); + PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py()); - // - // Loop over pad hits - for (AliRICHPadHit* mPad= - (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits); - mPad; - mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits)) + + + if (TMath::Abs(particle) < 10000000) { - Int_t cathode = mPad->fCathode; // cathode number - Int_t ipx = mPad->fPadX; // pad number on X - Int_t ipy = mPad->fPadY; // pad number on Y - Int_t iqpad = mPad->fQpad; // charge per pad - // - // - //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad); - - Float_t thex, they, thez; - segmentation=iChamber->GetSegmentationModel(cathode); - segmentation->GetPadC(ipx,ipy,thex,they,thez); - new((*pAddress)[countadr++]) TVector(2); - TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]); - trinfo(0)=(Float_t)track; - trinfo(1)=(Float_t)iqpad; - - digits[0]=ipx; - digits[1]=ipy; - digits[2]=iqpad; - - AliRICHTransientDigit* pdigit; - // build the list of fired pads and update the info - if (!pHitMap[nch]->TestHit(ipx, ipy)) { - list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter); - pHitMap[nch]->SetHit(ipx, ipy, counter); - counter++; - pdigit=(AliRICHTransientDigit*)list->At(list->GetLast()); - // list of tracks - TObjArray *trlist=(TObjArray*)pdigit->TrackList(); - trlist->Add(&trinfo); - } else { - pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy); - // update charge - (*pdigit).fSignal+=iqpad; - // update list of tracks - TObjArray* trlist=(TObjArray*)pdigit->TrackList(); - Int_t lastEntry=trlist->GetLast(); - TVector *ptrkP=(TVector*)trlist->At(lastEntry); - TVector &ptrk=*ptrkP; - Int_t lastTrack=Int_t(ptrk(0)); - Int_t lastCharge=Int_t(ptrk(1)); - if (lastTrack==track) { - lastCharge+=iqpad; - trlist->RemoveAt(lastEntry); - trinfo(0)=lastTrack; - trinfo(1)=lastCharge; - trlist->AddAt(&trinfo,lastEntry); - } else { - trlist->Add(&trinfo); + hitsTheta->Fill(theta,(float) 1); + if (R<5) + { + if (PTvertex>.5 && PTvertex<=1) + { + hitsTheta500MeV->Fill(theta,(float) 1); + } + if (PTvertex>1 && PTvertex<=2) + { + hitsTheta1GeV->Fill(theta,(float) 1); + } + if (PTvertex>2 && PTvertex<=3) + { + hitsTheta2GeV->Fill(theta,(float) 1); + } + if (PTvertex>3) + { + hitsTheta3GeV->Fill(theta,(float) 1); + } } - // check the track list - Int_t nptracks=trlist->GetEntriesFast(); - if (nptracks > 2) { - printf("Attention - tracks: %d (>2)\n",nptracks); - //printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy); - for (Int_t tr=0;trAt(tr); - TVector &pptrk=*pptrkP; - trk[tr]=Int_t(pptrk(0)); - chtrk[tr]=Int_t(pptrk(1)); - } - } // end if nptracks - } // end if pdigit - } //end loop over clusters - }// track type condition - } // hit loop - } // track loop - - // open the file with background - - if (addBackground ) { - ntracks =(Int_t)TrH1->GetEntries(); - //printf("background - icat,ntracks1 %d %d\n",icat,ntracks); - //printf("background - Start loop over tracks \n"); - // - // Loop over tracks - // - for (Int_t trak=0; trakClear(); - if (fClusters2) fClusters2->Clear(); - TrH1->GetEvent(trak); - // - // Loop over hits - AliRICHHit* mHit; - for(int j=0;jGetEntriesFast();++j) - { - mHit=(AliRICHHit*) (*fHits2)[j]; - Int_t nch = mHit->fChamber-1; // chamber number - if (nch >6) continue; - iChamber = &(pRICH->Chamber(nch)); - Int_t rmin = (Int_t)iChamber->RInner(); - Int_t rmax = (Int_t)iChamber->ROuter(); - // - // Loop over pad hits - for (AliRICHPadHit* mPad= - (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2); - mPad; - mPad=(AliRICHPadHit*)pRICH->NextPad(fClusters2)) - { - Int_t cathode = mPad->fCathode; // cathode number - Int_t ipx = mPad->fPadX; // pad number on X - Int_t ipy = mPad->fPadY; // pad number on Y - Int_t iqpad = mPad->fQpad; // charge per pad - - Float_t thex, they, thez; - segmentation=iChamber->GetSegmentationModel(cathode); - segmentation->GetPadC(ipx,ipy,thex,they,thez); - Float_t rpad=TMath::Sqrt(thex*thex+they*they); - if (rpad < rmin || iqpad ==0 || rpad > rmax) continue; - new((*pAddress)[countadr++]) TVector(2); - TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]); - trinfo(0)=-1; // tag background - trinfo(1)=-1; - digits[0]=ipx; - digits[1]=ipy; - digits[2]=iqpad; - if (trak <4 && nch==0) - printf("bgr - pHitMap[nch]->TestHit(ipx, ipy),trak %d %d\n", - pHitMap[nch]->TestHit(ipx, ipy),trak); - AliRICHTransientDigit* pdigit; - // build the list of fired pads and update the info - if (!pHitMap[nch]->TestHit(ipx, ipy)) { - list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter); - - pHitMap[nch]->SetHit(ipx, ipy, counter); - counter++; - printf("bgr new elem in list - counter %d\n",counter); - pdigit=(AliRICHTransientDigit*)list->At(list->GetLast()); - // list of tracks - TObjArray *trlist=(TObjArray*)pdigit->TrackList(); - trlist->Add(&trinfo); - } else { - pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy); - // update charge - (*pdigit).fSignal+=iqpad; - // update list of tracks - TObjArray* trlist=(TObjArray*)pdigit->TrackList(); - Int_t lastEntry=trlist->GetLast(); - TVector *ptrkP=(TVector*)trlist->At(lastEntry); - TVector &ptrk=*ptrkP; - Int_t lastTrack=Int_t(ptrk(0)); - if (lastTrack==-1) { - continue; - } else { - trlist->Add(&trinfo); - } - // check the track list - Int_t nptracks=trlist->GetEntriesFast(); - if (nptracks > 0) { - for (Int_t tr=0;trAt(tr); - TVector &pptrk=*pptrkP; - trk[tr]=Int_t(pptrk(0)); - chtrk[tr]=Int_t(pptrk(1)); + } + + //if (nch == 3) + //{ + + //printf("Particle type: %d\n",current->GetPdgCode()); + if (TMath::Abs(particle) < 50000051) + { + //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112) + if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050) + { + //gMC->Rndm(&random, 1); + if (random->Rndm() < .1) + production->Fill(current->Vz(),R,(float) 1); + if (TMath::Abs(particle) == 50000050) + //if (TMath::Abs(particle) > 50000000) + { + photons +=1; + if (R<5) + { + primaryphotons +=1; + if (current->Energy()>0.001) + highprimaryphotons +=1; + } + } + if (TMath::Abs(particle) == 2112) + { + neutron +=1; + if (current->Energy()>0.0001) + highneutrons +=1; + } } - } // end if nptracks - } // end if pdigit - } //end loop over clusters - } // hit loop - } // track loop - TTree *fAli=gAlice->TreeK(); - if (fAli) pFile =fAli->GetCurrentFile(); - pFile->cd(); - } // if Add - - Int_t tracks[10]; - Int_t charges[10]; - //cout<<"Start filling digits \n "<GetEntriesFast(); - //printf(" \n \n nentries %d \n",nentries); - - // start filling the digits - - for (Int_t nent=0;nentAt(nent); - if (address==0) continue; - - Int_t ich=address->fChamber; - Int_t q=address->fSignal; - iChamber=(AliRICHChamber*) (*fChambers)[ich]; - AliRICHResponse * response=iChamber->GetResponseModel(); - Int_t adcmax= (Int_t) response->MaxAdc(); - - - // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2) - //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY)); - Int_t pedestal = iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY); + if (TMath::Abs(particle) < 50000000) + { + production->Fill(current->Vz(),R,(float) 1); + //printf("Adding %d at %f\n",particle,R); + } + //mip->Fill(x,y,(float) 1); + } + + if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111) + { + if (R<5) + { + pionptspectravertex->Fill(PTvertex,(float) 1); + pionptspectrafinal->Fill(PTfinal,(float) 1); + } + } + + if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 + || TMath::Abs(particle)==311) + { + if (R<5) + { + kaonptspectravertex->Fill(PTvertex,(float) 1); + kaonptspectrafinal->Fill(PTfinal,(float) 1); + } + } + + + if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111) + { + pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf ("fParticle: %d, PDG code:%d\n",particle,current->GetPdgCode()); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + { + pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R); + } + //printf("Pion mass: %e\n",current->GetCalcMass()); + pion +=1; + if (TMath::Abs(particle)==211) + { + chargedpions +=1; + if (R<5) + { + primarypions +=1; + if (current->Energy()>1) + highprimarypions +=1; + } + } + } + if (TMath::Abs(particle)==2212) + { + protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //ptspectra->Fill(Pt,(float) 1); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass()); + proton +=1; + } + if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 + || TMath::Abs(particle)==311) + { + kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //ptspectra->Fill(Pt,(float) 1); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf("Kaon mass: %e\n",current->GetCalcMass()); + kaon +=1; + if (TMath::Abs(particle)==321) + { + chargedkaons +=1; + if (R<5) + { + primarykaons +=1; + if (current->Energy()>1) + highprimarykaons +=1; + } + } + } + if (TMath::Abs(particle)==11) + { + electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //ptspectra->Fill(Pt,(float) 1); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf("Electron mass: %e\n",current->GetCalcMass()); + if (particle == 11) + electron +=1; + if (particle == -11) + positron +=1; + } + if (TMath::Abs(particle)==13) + { + muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //ptspectra->Fill(Pt,(float) 1); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf("Muon mass: %e\n",current->GetCalcMass()); + muon +=1; + } + if (TMath::Abs(particle)==2112) + { + neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //ptspectra->Fill(Pt,(float) 1); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + { + neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R); + } + //printf("Neutron mass: %e\n",current->GetCalcMass()); + neutron +=1; + } + if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321) + { + if (current->Energy()-current->GetCalcMass()>1) + { + chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5) + chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + if (R>250 && R<450) + chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1); + } + } + //printf("Hits:%d\n",hit); + //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y); + // Fill the histograms + //Nh1+=nhits; + //h->Fill(x,y,(float) 1); + //} + //} + } + + } + + } + // } - //printf("Pedestal:%d\n",pedestal); - //Int_t pedestal=0; - Float_t treshold = (pedestal + 4*2.2); - - Float_t meanNoise = gRandom->Gaus(2.2, 0.3); - Float_t noise = gRandom->Gaus(0, meanNoise); - q+=(Int_t)(noise + pedestal); - //q+=(Int_t)(noise); - // magic number to be parametrised !!! - if ( q <= treshold) - { - q = q - pedestal; - continue; - } - q = q - pedestal; - if ( q >= adcmax) q=adcmax; - digits[0]=address->fPadX; - digits[1]=address->fPadY; - digits[2]=q; - - TObjArray* trlist=(TObjArray*)address->TrackList(); - Int_t nptracks=trlist->GetEntriesFast(); - - // this was changed to accomodate the real number of tracks - if (nptracks > 10) { - cout<<"Attention - tracks > 10 "< 2) { - printf("Attention - tracks > 2 %d \n",nptracks); - //printf("cat,ich,ix,iy,q %d %d %d %d %d \n", - //icat,ich,digits[0],digits[1],q); - } - for (Int_t tr=0;trAt(tr); - TVector &pp =*ppP; - tracks[tr]=Int_t(pp(0)); - charges[tr]=Int_t(pp(1)); - } //end loop over list of tracks for one pad - if (nptracks < 10 ) { - for (Int_t t=nptracks; t<10; t++) { - tracks[t]=0; - charges[t]=0; - } - } - //write file - if (ich==2) - fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]); - - // fill digits - pRICH->AddDigits(ich,tracks,charges,digits); - } - gAlice->TreeD()->Fill(); - - list->Delete(); - for(Int_t ii=0;iiSetPalette(1,0); + mystyle->cd(); + + //Create canvases, set the view range, show histograms + + TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150); + c2->Divide(2,2); + //c2->SetFillColor(42); + + c2->cd(1); + hitsTheta500MeV->SetFillColor(5); + hitsTheta500MeV->Draw(); + c2->cd(2); + hitsTheta1GeV->SetFillColor(5); + hitsTheta1GeV->Draw(); + c2->cd(3); + hitsTheta2GeV->SetFillColor(5); + hitsTheta2GeV->Draw(); + c2->cd(4); + hitsTheta3GeV->SetFillColor(5); + hitsTheta3GeV->Draw(); - //TTree *TD=gAlice->TreeD(); - //Stat_t ndig=TD->GetEntries(); - //cout<<"number of digits "<DigitsAddress(k); - int ndigit=fDch->GetEntriesFast(); - printf ("Chamber %d digits %d \n",k,ndigit); - } - pRICH->ResetDigits(); - char hname[30]; - sprintf(hname,"TreeD%d",nev); - gAlice->TreeD()->Write(hname,kOverwrite,0); - - // reset tree - // gAlice->TreeD()->Reset(); - delete list; - pAddress->Clear(); - // gObjectTable->Print(); + + + TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600); + c15->cd(); + production->SetFillColor(42); + production->SetXTitle("z (m)"); + production->SetYTitle("R (m)"); + production->Draw(); + + TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700); + c10->Divide(2,2); + c10->cd(1); + pionptspectravertex->SetFillColor(5); + pionptspectravertex->SetXTitle("Pt (GeV)"); + pionptspectravertex->Draw(); + c10->cd(2); + pionptspectrafinal->SetFillColor(5); + pionptspectrafinal->SetXTitle("Pt (GeV)"); + pionptspectrafinal->Draw(); + c10->cd(3); + kaonptspectravertex->SetFillColor(5); + kaonptspectravertex->SetXTitle("Pt (GeV)"); + kaonptspectravertex->Draw(); + c10->cd(4); + kaonptspectrafinal->SetFillColor(5); + kaonptspectrafinal->SetXTitle("Pt (GeV)"); + kaonptspectrafinal->Draw(); + + + TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350); + c16->Divide(2,1); + + c16->cd(1); + //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700); + electronspectra1->SetFillColor(5); + electronspectra1->SetXTitle("log(GeV)"); + electronspectra2->SetFillColor(46); + electronspectra2->SetXTitle("log(GeV)"); + electronspectra3->SetFillColor(10); + electronspectra3->SetXTitle("log(GeV)"); + //c13->SetLogx(); + electronspectra1->Draw(); + electronspectra2->Draw("same"); + electronspectra3->Draw("same"); + + c16->cd(2); + //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700); + muonspectra1->SetFillColor(5); + muonspectra1->SetXTitle("log(GeV)"); + muonspectra2->SetFillColor(46); + muonspectra2->SetXTitle("log(GeV)"); + muonspectra3->SetFillColor(10); + muonspectra3->SetXTitle("log(GeV)"); + //c14->SetLogx(); + muonspectra1->Draw(); + muonspectra2->Draw("same"); + muonspectra3->Draw("same"); + + //c16->cd(3); + //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700); + //neutronspectra1->SetFillColor(42); + //neutronspectra1->SetXTitle("log(GeV)"); + //neutronspectra2->SetFillColor(46); + //neutronspectra2->SetXTitle("log(GeV)"); + //neutronspectra3->SetFillColor(10); + //neutronspectra3->SetXTitle("log(GeV)"); + //c16->SetLogx(); + //neutronspectra1->Draw(); + //neutronspectra2->Draw("same"); + //neutronspectra3->Draw("same"); + + TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700); + //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700); + c9->Divide(2,2); + + c9->cd(1); + pionspectra1->SetFillColor(5); + pionspectra1->SetXTitle("log(GeV)"); + pionspectra2->SetFillColor(46); + pionspectra2->SetXTitle("log(GeV)"); + pionspectra3->SetFillColor(10); + pionspectra3->SetXTitle("log(GeV)"); + //c9->SetLogx(); + pionspectra1->Draw(); + pionspectra2->Draw("same"); + pionspectra3->Draw("same"); + + c9->cd(2); + //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700); + protonspectra1->SetFillColor(5); + protonspectra1->SetXTitle("log(GeV)"); + protonspectra2->SetFillColor(46); + protonspectra2->SetXTitle("log(GeV)"); + protonspectra3->SetFillColor(10); + protonspectra3->SetXTitle("log(GeV)"); + //c10->SetLogx(); + protonspectra1->Draw(); + protonspectra2->Draw("same"); + protonspectra3->Draw("same"); + + c9->cd(3); + //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700); + kaonspectra1->SetFillColor(5); + kaonspectra1->SetXTitle("log(GeV)"); + kaonspectra2->SetFillColor(46); + kaonspectra2->SetXTitle("log(GeV)"); + kaonspectra3->SetFillColor(10); + kaonspectra3->SetXTitle("log(GeV)"); + //c11->SetLogx(); + kaonspectra1->Draw(); + kaonspectra2->Draw("same"); + kaonspectra3->Draw("same"); + + c9->cd(4); + //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700); + chargedspectra1->SetFillColor(5); + chargedspectra1->SetXTitle("log(GeV)"); + chargedspectra2->SetFillColor(46); + chargedspectra2->SetXTitle("log(GeV)"); + chargedspectra3->SetFillColor(10); + chargedspectra3->SetXTitle("log(GeV)"); + //c12->SetLogx(); + chargedspectra1->Draw(); + chargedspectra2->Draw("same"); + chargedspectra3->Draw("same"); + + + + printf("*****************************************\n"); + printf("* Particle * Counts *\n"); + printf("*****************************************\n"); + + printf("* Pions: * %4d *\n",pion); + printf("* Charged Pions: * %4d *\n",chargedpions); + printf("* Primary Pions: * %4d *\n",primarypions); + printf("* Primary Pions (p>1GeV/c): * %4d *\n",highprimarypions); + printf("* Kaons: * %4d *\n",kaon); + printf("* Charged Kaons: * %4d *\n",chargedkaons); + printf("* Primary Kaons: * %4d *\n",primarykaons); + printf("* Primary Kaons (p>1GeV/c): * %4d *\n",highprimarykaons); + printf("* Muons: * %4d *\n",muon); + printf("* Electrons: * %4d *\n",electron); + printf("* Positrons: * %4d *\n",positron); + printf("* Protons: * %4d *\n",proton); + printf("* All Charged: * %4d *\n",(chargedpions+chargedkaons+muon+electron+positron+proton)); + printf("*****************************************\n"); + //printf("* Photons: * %3.1f *\n",photons); + //printf("* Primary Photons: * %3.1f *\n",primaryphotons); + //printf("* Primary Photons (p>1MeV/c):* %3.1f *\n",highprimaryphotons); + //printf("*****************************************\n"); + //printf("* Neutrons: * %3.1f *\n",neutron); + //printf("* Neutrons (p>100keV/c): * %3.1f *\n",highneutrons); + //printf("*****************************************\n"); + + if (gAlice->TreeD()) + { + gAlice->TreeD()->GetEvent(0); + + Float_t occ[7]; + Float_t sum=0; + Float_t mean=0; + printf("\n*****************************************\n"); + printf("* Chamber * Digits * Occupancy *\n"); + printf("*****************************************\n"); + + for (Int_t ich=0;ich<7;ich++) + { + TClonesArray *Digits = DigitsAddress(ich); // Raw clusters branch + Int_t ndigits = Digits->GetEntriesFast(); + occ[ich] = Float_t(ndigits)/(160*144); + sum += Float_t(ndigits)/(160*144); + printf("* %d * %d * %3.1f%% *\n",ich,ndigits,occ[ich]*100); + } + mean = sum/7; + printf("*****************************************\n"); + printf("* Mean occupancy * %3.1f%% *\n",mean*100); + printf("*****************************************\n"); + } + + printf("\nEnd of analysis\n"); + } -AliRICH& AliRICH::operator=(const AliRICH& rhs) +//_________________________________________________________________________________________________ + + +void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2) { -// Assignment operator - return *this; + +AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH"); + AliRICHSegmentationV0* segmentation; + AliRICHChamber* chamber; + + chamber = &(pRICH->Chamber(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 -} + //Int_t Pad[144][160]; + /*for (Int_t i=0;iGetEvent(nev); + + + //cout<<"nev "<TreeH(); + Stat_t ntracks = TH->GetEntries(); + + // Start loop on tracks in the hits containers + //Int_t Nc=0; + for (Int_t track=0; trackResetHits(); + TH->GetEvent(track); + Int_t nhits = pRICH->Hits()->GetEntriesFast(); + if (nhits) Nh+=nhits; + printf("Hits : %d\n",nhits); + for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1); + mHit; + mHit=(AliRICHHit*)pRICH->NextHit()) + { + 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->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); + + //printf("Particle:%9d\n",particle); + + TParticle *current = (TParticle*)gAlice->Particle(index); + //printf("Particle type: %d\n",sizeoff(Particles)); + + hitsTheta->Fill(theta,(float) 1); + //hitsPhi->Fill(phi,(float) 1); + //if (pRICH->GetDebugLevel() == -1) + //printf("Theta:%f, Phi:%f\n",theta,phi); + + //printf("Debug Level:%d\n",pRICH->GetDebugLevel()); + + if (current->GetPdgCode() < 10000000) + { + mip->Fill(x,y,(float) 1); + //printf("adding mip\n"); + //if (current->Energy() - current->GetCalcMass()>1 && freon==1) + //{ + hitsPhi->Fill(TMath::Abs(phi),(float) 1); + //hitsTheta->Fill(theta,(float) 1); + //printf("Theta:%f, Phi:%f\n",theta,phi); + //} + } + + if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111) + { + pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1); + } + if (TMath::Abs(particle)==2212) + { + protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1); + } + if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 + || TMath::Abs(particle)==311) + { + kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1); + } + if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321) + { + if (current->Energy() - current->GetCalcMass()>1) + chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1); + } + //printf("Hits:%d\n",hit); + //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y); + // Fill the histograms + Nh1+=nhits; + h->Fill(x,y,(float) 1); + //} + //} + } + + Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast(); + //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040) + //totalphotonsevent->Fill(ncerenkovs,(float) 1); + + if (ncerenkovs) { + printf("Cerenkovs : %d\n",ncerenkovs); + totalphotonsevent->Fill(ncerenkovs,(float) 1); + for (Int_t hit=0;hitCerenkovs()->UncheckedAt(hit); + Int_t nchamber = cHit->fChamber; // chamber number + Int_t index = cHit->Track(); + //Int_t pindex = (Int_t)(cHit->fIndex); + 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 + + 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); + Float_t energyckov = current->Energy(); + + if (current->GetPdgCode() == 50000051) + { + if (closs==4) + { + feedback->Fill(cx,cy,(float) 1); + feed++; + } + } + if (current->GetPdgCode() == 50000050) + { + + if (closs !=4) + { + phspectra2->Fill(energyckov*1e9,(float) 1); + } + + if (closs==4) + { + cerenkov->Fill(cx,cy,(float) 1); + + //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy); + + //TParticle *MIP = (TParticle*)gAlice->Particle(cmother); + AliRICHHit* mipHit = (AliRICHHit*) pRICH->Hits()->UncheckedAt(0); + mom[0] = current->Px(); + mom[1] = current->Py(); + mom[2] = current->Pz(); + //mom[0] = cHit->fMomX; + // mom[1] = cHit->fMomZ; + //mom[2] = cHit->fMomY; + //Float_t energymip = MIP->Energy(); + //Float_t Mip_px = mipHit->fMomFreoX; + //Float_t Mip_py = mipHit->fMomFreoY; + //Float_t Mip_pz = mipHit->fMomFreoZ; + //Float_t Mip_px = MIP->Px(); + //Float_t Mip_py = MIP->Py(); + //Float_t Mip_pz = MIP->Pz(); + + + + //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2]; + //Float_t rt = TMath::Sqrt(r); + //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz; + //Float_t Mip_rt = TMath::Sqrt(Mip_r); + //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001); + //Float_t cherenkov = TMath::ACos(coscerenkov); + ckovangle->Fill(cherenkov,(float) 1); //Cerenkov angle calculus + //printf("Cherenkov: %f\n",cherenkov); + Float_t ckphi=TMath::ATan2(mom[0], mom[2]); + hckphi->Fill(ckphi,(float) 1); + + + //Float_t mix = MIP->Vx(); + //Float_t miy = MIP->Vy(); + 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 = 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); + radius->Fill(final_radius,(float) 1); + + phspectra1->Fill(energyckov*1e9,(float) 1); + phot++; + } + for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){ + if (cmother == nmothers){ + if (closs == 4) + mothers2[cmother]++; + mothers[cmother]++; + } + } + } + } + } + + + if(gAlice->TreeR()) + { + Int_t nent=(Int_t)gAlice->TreeR()->GetEntries(); + gAlice->TreeR()->GetEvent(nent-1); + TClonesArray *Rawclusters = pRICH->RawClustAddress(2); // Raw clusters branch + //printf ("Rawclusters:%p",Rawclusters); + Int_t nrawclusters = Rawclusters->GetEntriesFast(); + + if (nrawclusters) { + printf("Raw Clusters : %d\n",nrawclusters); + for (Int_t hit=0;hitRawClustAddress(2)->UncheckedAt(hit); + //Int_t nchamber = rcHit->fChamber; // chamber number + //Int_t nhit = cHit->fHitNumber; // hit number + Int_t qtot = rcHit->fQ; // charge + Float_t fx = rcHit->fX; // x-position + Float_t fy = rcHit->fY; // y-position + //Int_t type = rcHit->fCtype; // cluster type ? + Int_t mult = rcHit->fMultiplicity; // How many pads form the cluster + pads += mult; + if (qtot > 0) { + //printf ("fx: %d, fy: %d\n",fx,fy); + if (fx>(x-4) && fx<(x+4) && fy>(y-4) && fy<(y+4)) { + //printf("There %d \n",mult); + padmip+=mult; + } else { + padnumber->Fill(mult,(float) 1); + nraw++; + if (mult<4) Clcharge->Fill(qtot,(float) 1); + } + + } + } + } + + + TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2); + Int_t nrechits1D = RecHits1D->GetEntriesFast(); + //printf (" nrechits:%d\n",nrechits); + + if(nrechits1D) + { + for (Int_t hit=0;hitRecHitsAddress1D(2)->UncheckedAt(hit); + Float_t r_omega = recHit1D->fOmega; // Cerenkov angle + Float_t *cer_pho = recHit1D->fCerPerPhoton; // Cerenkov angle per photon + Int_t *padsx = recHit1D->fPadsUsedX; // Pads Used fo reconstruction (x) + Int_t *padsy = recHit1D->fPadsUsedY; // Pads Used fo reconstruction (y) + Int_t goodPhotons = recHit1D->fGoodPhotons; // Number of pads used for reconstruction + + Omega1D->Fill(r_omega,(float) 1); + + for (Int_t i=0; iFill(cer_pho[i],(float) 1); + PadsUsed->Fill(padsx[i],padsy[i],1); + //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]); + } + + //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi); + } + } + + + TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2); + Int_t nrechits3D = RecHits3D->GetEntriesFast(); + //printf (" nrechits:%d\n",nrechits); + + if(nrechits3D) + { + 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]); + } + + clusev->Fill(nraw,(float) 1); + photev->Fill(phot,(float) 1); + feedev->Fill(feed,(float) 1); + padsmip->Fill(padmip,(float) 1); + padscl->Fill(pads,(float) 1); + //printf("Photons:%d\n",phot); + phot = 0; + feed = 0; + 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->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++) + { + TClonesArray *Digits = pRICH->DigitsAddress(ich); // Raw clusters branch + Int_t ndigits = Digits->GetEntriesFast(); + //printf("Digits:%d\n",ndigits); + padsev->Fill(ndigits,(float) 1); + if (ndigits) { + for (Int_t hit=0;hitUncheckedAt(hit); + //Int_t nchamber = dHit->GetChamber(); // chamber number + //Int_t nhit = dHit->fHitNumber; // hit number + 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); + if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot); + if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot); + } + } + } + } + } + + 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 + + TCanvas *c1 = 0; + TCanvas *c2 = 0; + TCanvas *c3 = 0; + TCanvas *c4 = 0; + TCanvas *c5 = 0; + TCanvas *c6 = 0; + TCanvas *c7 = 0; + TCanvas *c8 = 0; + TCanvas *c9 = 0; + 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(); + TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2); + Int_t nrechits1D = RecHits1D->GetEntriesFast(); + + switch(diaglevel) + { + case 1: + + c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350); + hc0->SetXTitle("ix (npads)"); + hc0->Draw("colz"); + // -// Calls the charge disintegration method of the current chamber and adds -// the simulated cluster to the root treee + c2 = new TCanvas("c2","Hits per type",100,100,600,700); + c2->Divide(2,2); + //c4->SetFillColor(42); + + c2->cd(1); + feedback->SetXTitle("x (cm)"); + feedback->SetYTitle("y (cm)"); + feedback->Draw("colz"); + + c2->cd(2); + //mip->SetFillColor(5); + mip->SetXTitle("x (cm)"); + mip->SetYTitle("y (cm)"); + mip->Draw("colz"); + + c2->cd(3); + //cerenkov->SetFillColor(5); + cerenkov->SetXTitle("x (cm)"); + cerenkov->SetYTitle("y (cm)"); + cerenkov->Draw("colz"); + + c2->cd(4); + //h->SetFillColor(5); + h->SetXTitle("x (cm)"); + h->SetYTitle("y (cm)"); + h->Draw("colz"); + + c3 = new TCanvas("c3","Hits distribution",150,150,600,350); + c3->Divide(2,1); + //c10->SetFillColor(42); + + c3->cd(1); + hitsX->SetFillColor(5); + hitsX->SetXTitle("(cm)"); + hitsX->Draw(); + + c3->cd(2); + hitsY->SetFillColor(5); + hitsY->SetXTitle("(cm)"); + hitsY->Draw(); + + + break; // - Int_t clhits[kNCH]; - Float_t newclust[6][500]; - Int_t nnew; - + case 2: + + c4 = new TCanvas("c4","Photon Spectra",50,50,600,350); + c4->Divide(2,1); + //c6->SetFillColor(42); + + c4->cd(1); + phspectra2->SetFillColor(5); + phspectra2->SetXTitle("energy (eV)"); + phspectra2->Draw(); + c4->cd(2); + phspectra1->SetFillColor(5); + phspectra1->SetXTitle("energy (eV)"); + phspectra1->Draw(); + + c5 = new TCanvas("c5","Particles Spectra",100,100,600,700); + c5->Divide(2,2); + //c9->SetFillColor(42); + + c5->cd(1); + pionspectra->SetFillColor(5); + pionspectra->SetXTitle("(GeV)"); + pionspectra->Draw(); + + c5->cd(2); + protonspectra->SetFillColor(5); + protonspectra->SetXTitle("(GeV)"); + protonspectra->Draw(); + + c5->cd(3); + kaonspectra->SetFillColor(5); + kaonspectra->SetXTitle("(GeV)"); + kaonspectra->Draw(); + + c5->cd(4); + chargedspectra->SetFillColor(5); + chargedspectra->SetXTitle("(GeV)"); + chargedspectra->Draw(); + + break; + + case 3: + + + if(gAlice->TreeR()) + { + c6=new TCanvas("c6","Clusters Statistics",50,50,600,700); + c6->Divide(2,2); + //c3->SetFillColor(42); + + c6->cd(1); + //TPad* c6_1; + //c6_1->SetLogy(); + Clcharge->SetFillColor(5); + Clcharge->SetXTitle("ADC counts"); + if (evNumber2>10) + { + Clcharge->Fit("expo"); + //expo->SetLineColor(2); + //expo->SetLineWidth(3); + } + Clcharge->Draw(); + + c6->cd(2); + padnumber->SetFillColor(5); + padnumber->SetXTitle("(counts)"); + padnumber->Draw(); + + c6->cd(3); + clusev->SetFillColor(5); + clusev->SetXTitle("(counts)"); + if (evNumber2>10) + { + clusev->Fit("gaus"); + //gaus->SetLineColor(2); + //gaus->SetLineWidth(3); + } + clusev->Draw(); + + c6->cd(4); + padsmip->SetFillColor(5); + padsmip->SetXTitle("(counts)"); + padsmip->Draw(); + } + + if(evNumber2<1) + { + c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700); + mother->SetFillColor(5); + mother->SetXTitle("counts"); + mother->Draw(); + } + + c7 = new TCanvas("c7","Production Statistics",100,100,600,700); + c7->Divide(2,2); + //c7->SetFillColor(42); + + c7->cd(1); + totalphotonsevent->SetFillColor(5); + totalphotonsevent->SetXTitle("Photons (counts)"); + if (evNumber2>10) + { + totalphotonsevent->Fit("gaus"); + //gaus->SetLineColor(2); + //gaus->SetLineWidth(3); + } + totalphotonsevent->Draw(); + + c7->cd(2); + photev->SetFillColor(5); + photev->SetXTitle("(counts)"); + if (evNumber2>10) + { + photev->Fit("gaus"); + //gaus->SetLineColor(2); + //gaus->SetLineWidth(3); + } + photev->Draw(); + + c7->cd(3); + feedev->SetFillColor(5); + feedev->SetXTitle("(counts)"); + if (evNumber2>10) + { + feedev->Fit("gaus"); + //gaus->SetLineColor(2); + //gaus->SetLineWidth(3); + } + feedev->Draw(); + + c7->cd(4); + padsev->SetFillColor(5); + padsev->SetXTitle("(counts)"); + if (evNumber2>10) + { + padsev->Fit("gaus"); + //gaus->SetLineColor(2); + //gaus->SetLineWidth(3); + } + padsev->Draw(); + + break; + + case 4: + + + if(nrechits3D) + { + 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); + OriginalPhi->SetFillColor(5); + if (evNumber2>10) + OriginalPhi->Fit("gaus"); + OriginalPhi->Draw(); + + //recontructed data + c8->cd(3); + Phi->SetFillColor(5); + if (evNumber2>10) + Phi->Fit("gaus"); + Phi->Draw(); + + 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(); + + 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(); + + + 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(); + + } + + if(nrechits1D) + { + c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700); + c9->Divide(3,2); + //c5->SetFillColor(42); + + c9->cd(1); + ckovangle->SetFillColor(5); + ckovangle->SetXTitle("angle (radians)"); + ckovangle->Draw(); + + c9->cd(2); + radius->SetFillColor(5); + radius->SetXTitle("radius (cm)"); + radius->Draw(); + + c9->cd(3); + hc0->SetXTitle("pads"); + hc0->Draw("box"); + + c9->cd(5); + Omega1D->SetFillColor(5); + Omega1D->SetXTitle("angle (radians)"); + Omega1D->Draw(); + + c9->cd(4); + PhotonCer->SetFillColor(5); + PhotonCer->SetXTitle("angle (radians)"); + PhotonCer->Draw(); + + c9->cd(6); + PadsUsed->SetXTitle("pads"); + PadsUsed->Draw("box"); + } + + break; + + case 5: + + printf("Drawing histograms.../n"); + + //if (ndigits) + //{ + c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700); + c1->Divide(4,2); + //c1->SetFillColor(42); + + c10->cd(1); + hc1->SetXTitle("ix (npads)"); + hc1->Draw("box"); + c10->cd(2); + hc2->SetXTitle("ix (npads)"); + hc2->Draw("box"); + c10->cd(3); + hc3->SetXTitle("ix (npads)"); + hc3->Draw("box"); + c10->cd(4); + hc4->SetXTitle("ix (npads)"); + hc4->Draw("box"); + c10->cd(5); + hc5->SetXTitle("ix (npads)"); + hc5->Draw("box"); + c10->cd(6); + hc6->SetXTitle("ix (npads)"); + hc6->Draw("box"); + c10->cd(7); + hc7->SetXTitle("ix (npads)"); + hc7->Draw("box"); + c10->cd(8); + hc0->SetXTitle("ix (npads)"); + hc0->Draw("box"); + //} // -// Integrated pulse height on chamber - - clhits[0]=fNhits+1; - - ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res); - Int_t ic=0; + c11 = new TCanvas("c11","Hits per type",100,100,600,700); + c11->Divide(2,2); + //c4->SetFillColor(42); + + c11->cd(1); + feedback->SetXTitle("x (cm)"); + feedback->SetYTitle("y (cm)"); + feedback->Draw(); + + c11->cd(2); + //mip->SetFillColor(5); + mip->SetXTitle("x (cm)"); + mip->SetYTitle("y (cm)"); + mip->Draw(); + + c11->cd(3); + //cerenkov->SetFillColor(5); + cerenkov->SetXTitle("x (cm)"); + cerenkov->SetYTitle("y (cm)"); + cerenkov->Draw(); + + c11->cd(4); + //h->SetFillColor(5); + h->SetXTitle("x (cm)"); + h->SetYTitle("y (cm)"); + h->Draw(); + + c12 = new TCanvas("c12","Hits distribution",150,150,600,350); + c12->Divide(2,1); + //c10->SetFillColor(42); + + c12->cd(1); + hitsX->SetFillColor(5); + hitsX->SetXTitle("(cm)"); + hitsX->Draw(); + + c12->cd(2); + hitsY->SetFillColor(5); + hitsY->SetXTitle("(cm)"); + hitsY->Draw(); + + break; + + } + + + // calculate the number of pads which give a signal + + + //Int_t Np=0; + /*for (Int_t i=0;i< NpadX;i++) { + for (Int_t j=0;j< NpadY;j++) { + if (Pad[i][j]>=6){ + Np+=1; + } + } + }*/ + //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1); + 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]; -// -// Add new clusters - for (Int_t i=0; i 0) { - ic++; -// Cathode plane - clhits[1] = Int_t(newclust[5][i]); -// Cluster Charge - clhits[2] = Int_t(newclust[0][i]); -// Pad: ix - clhits[3] = Int_t(newclust[1][i]); -// Pad: iy - clhits[4] = Int_t(newclust[2][i]); -// Pad: charge - clhits[5] = Int_t(newclust[3][i]); -// Pad: chamber sector - clhits[6] = Int_t(newclust[4][i]); - - AddPadHit(clhits); - } + // + // one branch for digits per chamber + // + for (Int_t i=0; i