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