82463a228f4204a396798d329ff7476f935c5f70
[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   
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   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
36   if (fPMDLoader == 0x0)
37     {
38       printf("Can not find PMDLoader\n");
39     }
40
41
42   fPMDLoader->LoadHits("READ");
43
44   // This reads the PMD Hits tree and assigns the right track number
45   // to a cell and stores in the summable digits tree
46   //
47
48   const Int_t kPi0 = 111;
49   const Int_t kGamma = 22;
50   Int_t npmd;
51   Int_t trackno;
52   Int_t smnumber;
53   Int_t trackpid;
54   Int_t mtrackno;
55   Int_t mtrackpid;
56
57   Float_t xPos, yPos, zPos;
58   Int_t xpad = -1, ypad = -1;
59   Float_t edep;
60   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
61
62   for (Int_t ievt = 0; ievt < nevt; ievt++)
63     {
64
65       printf("Event Number = %d\n",ievt);
66       Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
67       printf("Number of Particles = %d\n",nparticles);
68       fRunLoader->GetEvent(ievt);
69       // ------------------------------------------------------- //
70       // Pointer to specific detector hits.
71       // Get pointers to Alice detectors and Hits containers
72       
73       TTree* treeH = fPMDLoader->TreeH();
74       
75       Int_t ntracks    = (Int_t) treeH->GetEntries();
76       printf("Number of Tracks in the TreeH = %d\n", ntracks);
77       
78       
79       TClonesArray* hits = 0;
80       if (fPMD) hits = fPMD->Hits();
81       
82       // Start loop on tracks in the hits containers
83       
84       for (Int_t track=0; track<ntracks;track++) 
85         {
86           gAlice->ResetHits();
87           treeH->GetEvent(track);
88           if (fPMD) 
89             {
90               npmd = hits->GetEntriesFast();
91               for (int ipmd = 0; ipmd < npmd; ipmd++) 
92                 {
93                   fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
94                   trackno = fPMDHit->GetTrack();
95
96                   //fprintf(fpw,"trackno = %d\n",trackno);
97
98                   //  get kinematics of the particles
99                   
100                   TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
101                   trackpid  = mparticle->GetPdgCode();
102                   
103                   Int_t igatr = -999;
104                   Int_t ichtr = -999;
105                   Int_t igapid = -999;
106                   Int_t imo;
107                   Int_t igen = 0;
108                   Int_t idmo = -999;
109                   
110                   Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
111                   if (mparticle->GetFirstMother() == -1)
112                     {
113                       tracknoOld  = trackno;
114                       trackpidOld = trackpid;
115                       statusOld   = -1;
116                       
117                       vx = mparticle->Vx();
118                       vy = mparticle->Vy();
119                       vz = mparticle->Vz();
120
121                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, -1, trackpid, vx, vy, vz);
122                       
123                     }
124                   Int_t igstatus = 0;
125                   while((imo = mparticle->GetFirstMother()) >= 0)
126                     {
127                       igen++;
128                       
129                       mparticle =  gAlice->GetMCApp()->Particle(imo);
130                       idmo = mparticle->GetPdgCode();
131                       
132                       vx = mparticle->Vx();
133                       vy = mparticle->Vy();
134                       vz = mparticle->Vz();
135                       
136                       //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
137                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
138                       
139                       if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
140                         {
141                           igatr = imo;
142                           igapid = idmo;
143                           igstatus = 1;
144                         }
145                       if(igstatus == 0)
146                         {
147                           if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
148                             {
149                               igatr = imo;
150                               igapid = idmo;
151                             }
152                         }
153                       ichtr = imo;
154                     } // end of while loop
155                   
156                   if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
157                     {
158                       mtrackno = igatr;
159                       mtrackpid = igapid;
160                     }
161                   else
162                     {
163                       mtrackno  = ichtr;
164                       mtrackpid = idmo;
165                     }
166                   if (statusOld == -1)
167                     {
168                       mtrackno  = tracknoOld;
169                       mtrackpid = trackpidOld;
170                     }
171                   
172                   xPos = fPMDHit->X();
173                   yPos = fPMDHit->Y();
174                   zPos = fPMDHit->Z();
175                   
176                   edep       = fPMDHit->GetEnergy();
177                   Int_t vol1 = fPMDHit->GetVolume(1); // Column
178                   Int_t vol2 = fPMDHit->GetVolume(2); // Row
179                   Int_t vol3 = fPMDHit->GetVolume(3); // UnitModule
180                   Int_t vol6 = fPMDHit->GetVolume(6); // SuperModule
181                   
182                   // -----------------------------------------//
183                   // For Super Module 1 & 2                   //
184                   //  nrow = 96, ncol = 48                    //
185                   // For Super Module 3 & 4                   //
186                   //  nrow = 48, ncol = 96                    //
187                   // -----------------------------------------//
188                   
189                   smnumber = (vol6-1)*6 + vol3;
190                   
191                   if (vol6 == 1 || vol6 == 2)
192                     {
193                       xpad = vol2;
194                       ypad = vol1;
195                     }
196                   else if (vol6 == 3 || vol6 == 4)
197                     {
198                       xpad = vol1;
199                       ypad = vol2;
200                     }
201                   
202                   //printf("Zposition = %f Edeposition = %f",zPos,edep);
203                   Float_t zposition = TMath::Abs(zPos);
204                   
205                 }
206             }
207         } // Track Loop ended
208       
209     }
210   fRunLoader->UnloadgAlice();
211   fRunLoader->UnloadHeader();
212   fRunLoader->UnloadKinematics();
213   
214   fPMDLoader->UnloadHits();
215
216   timer.Stop();
217   timer.Print();
218
219 }