fix
[u/mrichter/AliRoot.git] / PMD / AliPMDHitsRead.C
CommitLineData
515154c7 1//
2// This macro reads the Hits Tree
3//
4void AliPMDHitsRead(Int_t nevt = 1)
5{
6
7 TStopwatch timer;
8 timer.Start();
9ebfe36a 9 TH2F *h2 = new TH2F("h2"," Y vs. X",200,-100.,100.,200,-100.,100.);
515154c7 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");
585c05b4 35
515154c7 36 fPMDLoader = fRunLoader->GetLoader("PMDLoader");
585c05b4 37
515154c7 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;
9ebfe36a 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;
515154c7 59 Float_t edep;
60 Float_t vx = -999.0, vy = -999.0, vz = -999.0;
9ebfe36a 61 Float_t xPos, yPos, zPos;
62 Float_t xx, yy;
63
64 AliPMDUtility cc;
515154c7 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 {
585c05b4 90 gAlice->GetMCApp()->ResetHits();
515154c7 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
585c05b4 103
515154c7 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
585c05b4 176 //printf("mtrackno = %d mtrackpid = %d\n",mtrackno,mtrackpid);
177
515154c7 178 xPos = fPMDHit->X();
179 yPos = fPMDHit->Y();
180 zPos = fPMDHit->Z();
585c05b4 181
182 Float_t time = fPMDHit->GetTime();
515154c7 183
585c05b4 184 printf("++++++++++ time = %f\n",time);
185
515154c7 186 edep = fPMDHit->GetEnergy();
187 Int_t vol1 = fPMDHit->GetVolume(1); // Column
188 Int_t vol2 = fPMDHit->GetVolume(2); // Row
585c05b4 189
190 Int_t vol3 = fPMDHit->GetVolume(4); // UnitModule
191
515154c7 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
585c05b4 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
9ebfe36a 212 if(zPos > 361.5)
515154c7 213 {
9ebfe36a 214 cc.RectGeomCellPos(smnumber,xpad,ypad,xx,yy);
215 h2->Fill(xx,yy);
515154c7 216 }
217
515154c7 218 }
219 }
220 } // Track Loop ended
221
222 }
9ebfe36a 223
224 h2->Draw();
225
515154c7 226 fRunLoader->UnloadgAlice();
227 fRunLoader->UnloadHeader();
228 fRunLoader->UnloadKinematics();
515154c7 229 fPMDLoader->UnloadHits();
230
231 timer.Stop();
232 timer.Print();
233
234}