da documentation page is defined in the Link
[u/mrichter/AliRoot.git] / PMD / AliPMDHitsRead.C
1 //
2 // This macro reads the Hits Tree
3 //
4 void AliPMDHitsRead(Int_t nevt = 1)
5 {
6   
7   TStopwatch timer;
8   timer.Start();
9   TH2F *h2 = new TH2F("h2"," Y vs. X",200,-100.,100.,200,-100.,100.); 
10   //  FILE *fpw = fopen("alipmdhits.dat","w");
11
12   AliRunLoader *fRunLoader = AliRunLoader::Open("galice.root");
13
14   if (!fRunLoader)
15    {
16      printf("Can not open session for file ");
17    }
18   
19   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
20   if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
21   if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
22
23   gAlice = fRunLoader->GetAliRun();
24   
25   if (gAlice)
26     {
27       printf("Alirun object found\n");
28     }
29   else
30     {
31       printf("Could not found Alirun object\n");
32     }
33   
34   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
35
36   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
37
38   if (fPMDLoader == 0x0)
39     {
40       printf("Can not find PMDLoader\n");
41     }
42
43
44   fPMDLoader->LoadHits("READ");
45
46   // This reads the PMD Hits tree and assigns the right track number
47   // to a cell and stores in the summable digits tree
48   //
49
50   const Int_t kPi0 = 111;
51   const Int_t kGamma = 22;
52   Int_t   npmd;
53   Int_t   trackno;
54   Int_t   smnumber;
55   Int_t   trackpid;
56   Int_t   mtrackno;
57   Int_t   mtrackpid;
58   Int_t   xpad = -1, ypad = -1;
59   Float_t edep;
60   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
61   Float_t xPos, yPos, zPos;
62   Float_t xx, yy;
63
64   AliPMDUtility cc;
65
66   for (Int_t ievt = 0; ievt < nevt; ievt++)
67     {
68
69       printf("Event Number = %d\n",ievt);
70       Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
71       printf("Number of Particles = %d\n",nparticles);
72       fRunLoader->GetEvent(ievt);
73       // ------------------------------------------------------- //
74       // Pointer to specific detector hits.
75       // Get pointers to Alice detectors and Hits containers
76       
77       TTree* treeH = fPMDLoader->TreeH();
78       
79       Int_t ntracks    = (Int_t) treeH->GetEntries();
80       printf("Number of Tracks in the TreeH = %d\n", ntracks);
81       
82       
83       TClonesArray* hits = 0;
84       if (fPMD) hits = fPMD->Hits();
85       
86       // Start loop on tracks in the hits containers
87       
88       for (Int_t track=0; track<ntracks;track++) 
89         {
90           gAlice->GetMCApp()->ResetHits();
91           treeH->GetEvent(track);
92           if (fPMD) 
93             {
94               npmd = hits->GetEntriesFast();
95               for (int ipmd = 0; ipmd < npmd; ipmd++) 
96                 {
97                   fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
98                   trackno = fPMDHit->GetTrack();
99
100                   //fprintf(fpw,"trackno = %d\n",trackno);
101
102                   //  get kinematics of the particles
103
104                   TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
105                   trackpid  = mparticle->GetPdgCode();
106                   
107                   Int_t igatr = -999;
108                   Int_t ichtr = -999;
109                   Int_t igapid = -999;
110                   Int_t imo;
111                   Int_t igen = 0;
112                   Int_t idmo = -999;
113                   
114                   Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
115                   if (mparticle->GetFirstMother() == -1)
116                     {
117                       tracknoOld  = trackno;
118                       trackpidOld = trackpid;
119                       statusOld   = -1;
120                       
121                       vx = mparticle->Vx();
122                       vy = mparticle->Vy();
123                       vz = mparticle->Vz();
124
125                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, -1, trackpid, vx, vy, vz);
126                       
127                     }
128                   Int_t igstatus = 0;
129                   while((imo = mparticle->GetFirstMother()) >= 0)
130                     {
131                       igen++;
132                       
133                       mparticle =  gAlice->GetMCApp()->Particle(imo);
134                       idmo = mparticle->GetPdgCode();
135                       
136                       vx = mparticle->Vx();
137                       vy = mparticle->Vy();
138                       vz = mparticle->Vz();
139                       
140                       //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
141                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
142                       
143                       if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
144                         {
145                           igatr = imo;
146                           igapid = idmo;
147                           igstatus = 1;
148                         }
149                       if(igstatus == 0)
150                         {
151                           if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
152                             {
153                               igatr = imo;
154                               igapid = idmo;
155                             }
156                         }
157                       ichtr = imo;
158                     } // end of while loop
159                   
160                   if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
161                     {
162                       mtrackno = igatr;
163                       mtrackpid = igapid;
164                     }
165                   else
166                     {
167                       mtrackno  = ichtr;
168                       mtrackpid = idmo;
169                     }
170                   if (statusOld == -1)
171                     {
172                       mtrackno  = tracknoOld;
173                       mtrackpid = trackpidOld;
174                     }
175                   
176                   //printf("mtrackno = %d  mtrackpid = %d\n",mtrackno,mtrackpid);
177
178                   xPos = fPMDHit->X();
179                   yPos = fPMDHit->Y();
180                   zPos = fPMDHit->Z();
181
182                   Float_t time = fPMDHit->GetTime();
183                   
184                   printf("++++++++++ time = %f\n",time);
185
186                   edep       = fPMDHit->GetEnergy();
187                   Int_t vol1 = fPMDHit->GetVolume(1); // Column
188                   Int_t vol2 = fPMDHit->GetVolume(2); // Row
189
190                   Int_t vol3 = fPMDHit->GetVolume(4); // UnitModule
191
192                   
193                   // -----------------------------------------//
194                   // For Super Module 1 & 2                   //
195                   //  nrow = 96, ncol = 48                    //
196                   // For Super Module 3 & 4                   //
197                   //  nrow = 48, ncol = 96                    //
198                   // -----------------------------------------//
199                   
200                   if (vol3 < 24)
201                     {
202                       smnumber = vol3;
203                     }
204                   else
205                     {
206                       smnumber = vol3 - 24;
207                     }
208
209                   xpad = vol1;
210                   ypad = vol2;
211
212                   if(zPos > 361.5)
213                     {
214                       cc.RectGeomCellPos(smnumber,xpad,ypad,xx,yy);
215                       h2->Fill(xx,yy);
216                     }
217                   
218                 }
219             }
220         } // Track Loop ended
221       
222     }
223
224   h2->Draw();
225
226   fRunLoader->UnloadgAlice();
227   fRunLoader->UnloadHeader();
228   fRunLoader->UnloadKinematics();
229   fPMDLoader->UnloadHits();
230
231   timer.Stop();
232   timer.Print();
233
234 }