]>
Commit | Line | Data |
---|---|---|
23efe5f1 | 1 | |
2 | Int_t max_modules=2269; | |
3 | ||
4 | void dEdXyy (Int_t evNumber1=0,Int_t evNumber2=0,AliITSPid* pid=0,int mtest=0) | |
5 | { | |
6 | if (gClassTable->GetID("AliRun") < 0) { | |
7 | gROOT->LoadMacro("loadlibs.C"); | |
8 | loadlibs(); | |
9 | } else { | |
10 | delete gAlice; | |
11 | gAlice=0; | |
12 | } | |
13 | ||
14 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
15 | ||
16 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
17 | printf("file %p\n",file); | |
18 | if (file) file->Close(); | |
19 | file = new TFile("galice.root","UPDATE"); | |
20 | file->ls(); | |
21 | ||
22 | printf ("I'm after Map \n"); | |
23 | ||
24 | // Get AliRun object from file or create it if not on file | |
25 | ||
26 | if (!gAlice) { | |
27 | gAlice = (AliRun*)file->Get("gAlice"); | |
28 | if (gAlice) printf("AliRun object found on file\n"); | |
29 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
30 | } | |
31 | printf ("I'm after gAlice \n"); | |
32 | ||
33 | // Create Histogramms ------------------------------------------------- | |
34 | ||
35 | TH1F *dedxSDD = new TH1F("dedxSDD","Particle energy (KeV) for layer 3 and 4",100,0.,400.); | |
36 | //TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,400.); | |
37 | TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,3000.); | |
38 | //TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,6.); | |
39 | TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,30.); | |
40 | //TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,2000.); | |
41 | TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,15000.); | |
42 | //TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,6.); | |
43 | TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,30.); | |
44 | TH1F *dedxSSD = new TH1F("dedxSSD","Particle energy (KeV) for layer 5 and 6",100,0.,400.); | |
45 | //TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,400.); | |
46 | TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,3000.); | |
47 | //TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,6.); | |
48 | TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,30.); | |
49 | //TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,200.); | |
50 | TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,2000.); | |
51 | //TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,6.); | |
52 | TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,30.); | |
53 | ||
54 | ||
55 | ||
56 | /* | |
57 | TH1F *dedxSDD = new TH1F("dedxSDD","Particle energy (KeV) for layer 3 and 4",100,0.,4000.); | |
58 | TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,4000.); | |
59 | TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,40.); | |
60 | TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,20000.); | |
61 | TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,40.); | |
62 | TH1F *dedxSSD = new TH1F("dedxSSD","Particle energy (KeV) for layer 5 and 6",100,0.,4000.); | |
63 | TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,4000.); | |
64 | TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,40.); | |
65 | TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,2000.); | |
66 | TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,40.); | |
67 | */ | |
68 | ||
69 | TH1F *pathSDD = new TH1F("pathSDD","Path length in the SDD",100,0.,1000.); | |
70 | TH1F *pathSSD = new TH1F("pathSSD","Path length in the SSD",100,0.,1000.); | |
71 | ||
72 | // ------------------------------------------------------------- | |
73 | AliITS *ITS = (AliITS*) gAlice->GetModule("ITS"); | |
74 | if (!ITS) return; | |
75 | AliITSgeom *geom = ITS->GetITSgeom(); | |
76 | ||
77 | // Set the models for cluster finding | |
78 | ||
79 | // SPD | |
80 | ||
81 | ITS->MakeTreeC(); | |
82 | Int_t nparticles=gAlice->GetEvent(0); | |
83 | ||
84 | AliITSDetType *iDetType=ITS->DetType(0); | |
85 | AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); | |
86 | TClonesArray *dig0 = ITS->DigitsAddress(0); | |
87 | TClonesArray *recp0 = ITS->ClustersAddress(0); | |
88 | AliITSClusterFinderSPD *rec0=new AliITSClusterFinderSPD(seg0,dig0,recp0); | |
89 | ITS->SetReconstructionModel(0,rec0); | |
90 | ||
91 | // SDD | |
92 | ||
93 | AliITSDetType *iDetType=ITS->DetType(1); | |
94 | AliITSgeom *geom = ITS->GetITSgeom(); | |
95 | ||
96 | AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel(); | |
97 | if (!seg1) seg1 = new AliITSsegmentationSDD(geom); | |
98 | AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel(); | |
99 | if (!res1) res1=new AliITSresponseSDD(); | |
100 | ||
101 | Float_t baseline,noise; | |
102 | res1->GetNoiseParam(noise,baseline); | |
103 | Float_t noise_after_el = res1->GetNoiseAfterElectronics(); | |
104 | Float_t thres = baseline; | |
105 | thres += (4.*noise_after_el); // TB // (4.*noise_after_el); | |
106 | printf("thres %f\n",thres); | |
107 | //res1->Print(); | |
108 | ||
109 | TClonesArray *dig1 = ITS->DigitsAddress(1); | |
110 | TClonesArray *recp1 = ITS->ClustersAddress(1); | |
111 | AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1); | |
112 | rec1->SetCutAmplitude((int)thres); | |
113 | ITS->SetReconstructionModel(1,rec1); | |
114 | //rec1->Print(); | |
115 | ||
116 | // SSD | |
117 | ||
118 | AliITSDetType *iDetType=ITS->DetType(2); | |
119 | AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel(); | |
120 | TClonesArray *dig2 = ITS->DigitsAddress(2); | |
121 | AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2); | |
122 | ITS->SetReconstructionModel(2,rec2); | |
123 | ||
124 | // | |
125 | // Loop over events | |
126 | // | |
127 | Int_t Nh=0; | |
128 | Int_t Nh1=0; | |
129 | for (int nev=0; nev<= evNumber2; nev++) { | |
130 | Int_t nparticles = 0; | |
131 | nparticles = gAlice->GetEvent(nev); | |
132 | cout << "nev " << nev <<endl; | |
133 | cout << "nparticles " << nparticles <<endl; | |
134 | if (nev < evNumber1) continue; | |
135 | if (nparticles <= 0) return; | |
136 | ||
137 | AliITShit *itsHit; | |
138 | AliITShit *itsHitPrev; | |
139 | AliITSRecPoint *itsPnt = 0; | |
140 | AliITSRawClusterSDD *itsClu = 0; | |
141 | ||
142 | // Get Hit, Cluster & Recpoints Tree Pointers | |
143 | ||
144 | TTree *TH = gAlice->TreeH(); | |
145 | Int_t nenthit=TH->GetEntries(); | |
146 | printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit); | |
147 | ||
148 | ITS->GetTreeC(nev); | |
149 | TTree *TC=ITS->TreeC(); | |
150 | Int_t nentclu=TC->GetEntries(); | |
151 | printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu); | |
152 | ||
153 | TTree *TR = gAlice->TreeR(); | |
154 | Int_t nentrec=TR->GetEntries(); | |
155 | printf("Found %d entries in the RecPoints tree\n",nentrec); | |
156 | ||
157 | // Get Pointers to Clusters & Recpoints TClonesArrays | |
158 | ||
159 | TClonesArray *ITSclu = ITS->ClustersAddress(1); | |
160 | printf ("ITSclu %p \n",ITSclu); | |
161 | // Int_t ncl = ITSclu->GetEntries(); | |
162 | // cout<<"ncluster ="<<ncl<<endl; | |
163 | TClonesArray *ITSrec = ITS->RecPoints(); | |
164 | printf ("ITSrec %p \n",ITSrec); | |
165 | ||
166 | // check recpoints | |
167 | ||
168 | //Int_t nbytes = 0; | |
169 | Int_t totpoints = 0; | |
170 | Int_t totclust = 0; | |
171 | ||
172 | // check hits | |
173 | ||
174 | Int_t nmodules=0; | |
175 | ||
176 | ITS->InitModules(-1,nmodules); | |
177 | ITS->FillModules(nev,0,nmodules,"",""); | |
178 | ||
179 | TObjArray *fITSmodules = ITS->GetModules(); | |
180 | ||
181 | Int_t first0 = geom->GetStartDet(0); // SPD | |
182 | Int_t last0 = geom->GetLastDet(0); // SPD | |
183 | Int_t first1 = geom->GetStartDet(1); // SDD | |
184 | Int_t last1 = geom->GetLastDet(1); // SDD | |
185 | Int_t first2 = geom->GetStartDet(2); // SSD | |
186 | Int_t last2 = geom->GetLastDet(2); // SSD | |
187 | ||
188 | // For the SPD: first0 = 0, last0 = 239 (240 modules); | |
189 | // for the SDD: first1 = 240, last1 = 499 (260 modules); | |
190 | // for the SSD: first2 = 500, last2 = 2269 (1770 modules). | |
191 | ||
192 | Int_t mod; | |
193 | printf("det type %d first0, last0 %d %d \n",0,first0,last0); | |
194 | printf("det type %d first1, last1 %d %d \n",1,first1,last1); | |
195 | printf("det type %d first2, last2 %d %d \n",2,first2,last2); | |
196 | Int_t negtrSDD = 0; | |
197 | Int_t negtrSSD = 0; | |
198 | ||
199 | for (mod=0; mod<last2+1; mod++) { | |
200 | if(mod>max_modules)continue; | |
201 | if(mod < first1) continue; // for the SDD/SSD only | |
202 | ||
203 | AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod); | |
204 | Int_t nhits = Mod->GetNhits(); | |
205 | Int_t layer = Mod->GetLayer(); | |
206 | Int_t ladder = Mod->GetLadder(); | |
207 | Int_t detector = Mod->GetDet(); | |
208 | Int_t hit; | |
209 | Int_t parent; | |
210 | Int_t dray; | |
211 | Int_t hitstat; | |
212 | Int_t partcode; | |
213 | Int_t hitprim; | |
214 | Float_t epart; | |
215 | Float_t pathInSi=300; | |
216 | Float_t xhit; | |
217 | Float_t yhit; | |
218 | Float_t zhit; | |
219 | Float_t px; | |
220 | Float_t py; | |
221 | Float_t pz; | |
222 | Float_t pmod; | |
223 | Float_t xhit0 = 1e+7; | |
224 | Float_t yhit0 = 1e+7; | |
225 | Float_t zhit0 = 1e+7; | |
226 | ||
227 | TTree *TR = gAlice->TreeR(); | |
228 | ITS->ResetClusters(); | |
229 | //cout << "CLUSTERS: get" << endl; | |
230 | TC->GetEvent(mod); | |
231 | //cout << "RECPOINTS: reset" << endl; | |
232 | ITS->ResetRecPoints(); | |
233 | //cout << "RECPOINTS: get" << endl; | |
234 | //TR->GetEvent(mod+1); | |
235 | TR->GetEvent(mod); | |
236 | ||
237 | Int_t nrecp = ITSrec->GetEntries(); | |
238 | totpoints += nrecp; | |
239 | if (!nrecp) continue; | |
240 | ||
241 | Int_t trackRecp[3]; | |
242 | Int_t itr; | |
243 | Int_t TrackRecPoint; | |
244 | Float_t PmodRecP; | |
245 | Float_t signalRP; | |
246 | ||
247 | // cout <<" module,layer,ladder,detector,nrecp,nhits ="<< | |
248 | // mod<<","<<layer<<","<<ladder<<","<<detector<<","<<nrecp<<","<<nhits<< endl; | |
249 | cout<<"."; | |
250 | ||
251 | // ---------- RecPoint signal analysis for the SDD/SSD --------- | |
252 | ||
253 | for (Int_t pnt=0;pnt<nrecp;pnt++) { // recpoint loop | |
254 | ||
255 | itsPnt = (AliITSRecPoint*)ITSrec->At(pnt); | |
256 | ||
257 | Int_t RecPointPrim = 0; | |
258 | if(!itsPnt) continue; | |
259 | ||
260 | signalRP = itsPnt->GetQ(); | |
261 | trackRecp[0] = itsPnt->GetLabel(0); | |
262 | trackRecp[1] = itsPnt->GetLabel(1); | |
263 | trackRecp[2] = itsPnt->GetLabel(2); | |
264 | ||
265 | // cout<<"New Recp: pnt,signal,tr0,tr1,tr2 ="<<pnt | |
266 | // <<","<<signalRP<<","<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<endl; | |
267 | ||
268 | /* | |
269 | if(trackRecp[0]<0&&trackRecp[1]<0&&trackRecp[2]<0) { | |
270 | cout<<"pnt,tr0,1,2 ="<<pnt<<","<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<endl; | |
271 | if(mod<first2) negtrSDD += 1; | |
272 | if(mod>=first2) negtrSSD += 1; | |
273 | } | |
274 | */ | |
275 | ||
276 | xhit0 = 1e+7; | |
277 | yhit0 = 1e+7; | |
278 | zhit0 = 1e+7; | |
279 | ||
280 | // Hit loop | |
281 | for (hit=0;hit<nhits;hit++) { | |
282 | ||
283 | itsHit = (AliITShit*)Mod->GetHit(hit); | |
284 | ||
285 | hitprim = 0; | |
286 | dray = 0; | |
287 | Int_t hitRecp = 0; | |
288 | Int_t trackHit = itsHit->GetTrack(); | |
289 | hitstat = itsHit->GetTrackStatus(); | |
290 | zhit = 10000*itsHit->GetZL(); | |
291 | xhit = 10000*itsHit->GetXL(); | |
292 | yhit = 10000*itsHit->GetYL(); | |
293 | Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV | |
294 | ||
295 | // cout<<" New hit: hit,track,hitstat,yhit,ehit ="<<hit<<","<<trackHit<<","<<hitstat<<","<<yhit<<","<<dEn<<endl; | |
296 | ||
297 | // check the primary particle | |
298 | partcode = itsHit->GetParticle()->GetPdgCode(); | |
299 | parent = itsHit->GetParticle()->GetFirstMother(); | |
300 | if(parent < 0) hitprim = 1; // primery particle | |
301 | ||
302 | ||
303 | ||
304 | // select the primery hits with the track number equaled to | |
305 | // the one of the RecPoint | |
306 | for(itr=0;itr<3;itr++) { | |
307 | if(trackRecp[itr] == trackHit) hitRecp = 1; | |
308 | if(trackRecp[itr] == trackHit && hitprim == 1) { | |
309 | RecPointPrim = 1; | |
310 | TrackRecPoint = trackRecp[itr]; | |
311 | } | |
312 | if(trackRecp[itr] == trackHit && hitprim == 0) trackRecp[itr]=-100; | |
313 | } | |
314 | ||
315 | if(hitRecp != 1) continue; // hit doesn't correspond to the recpoint | |
316 | ||
317 | //cout<<"Hit Ok: trackRecp0,1,2,hitRecp,RecPointPrim,TrackRecPoint ="<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<","<<hitRecp<<","<<RecPointPrim<<","<<TrackRecPoint<<endl; | |
318 | ||
319 | // if(hitstat == 66 && yhit < -136.) { // entering hit | |
320 | if(hitstat == 66 && hitprim == 1) { // entering hit | |
321 | xhit0 = xhit; | |
322 | yhit0 = yhit; | |
323 | zhit0 = zhit; | |
324 | } | |
325 | // cout<<" xhit0,zhit0 ="<<xhit0<<","<<zhit0<<","<<endl; | |
326 | ||
327 | if(hitstat == 66) continue; // Take only the hits with the not | |
328 | // zero energy | |
329 | ||
330 | if(xhit0 > 9e+6 || zhit0 > 9e+6) { | |
331 | // cout<<"default xhit0,zhit0 ="<<xhit0<<","<<zhit0<<","<<endl; | |
332 | continue; | |
333 | } | |
334 | ||
335 | // check the particle code (type) | |
336 | // Int_t parent = itsHit->GetParticle()->GetFirstMother(); | |
337 | //partcode = itsHit->GetParticle()->GetPdgCode(); | |
338 | ||
339 | // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+ | |
340 | // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda | |
341 | ||
342 | /* | |
343 | px = itsHit->GetPXL(); // momentum for the hit | |
344 | py = itsHit->GetPYL(); | |
345 | pz = itsHit->GetPZL(); | |
346 | pmod = 1000*sqrt(px*px+py*py+pz*pz); | |
347 | */ | |
348 | ||
349 | pmod = itsHit->GetParticle()->P(); // momentum at the vertex | |
350 | pmod *= 1.0e+3; | |
351 | ||
352 | // cout<<"track,partcode,pmod,hitprim ="<<trackHit<<","<<partcode<<","<<pmod<<","<<hitprim<<endl; | |
353 | ||
354 | if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- | |
355 | // at p < 6 MeV/c | |
356 | ||
357 | // find the primery particle path length in Si and pmod (MeV/c) | |
358 | if((hitstat == 68 || hitstat == 33) && hitprim == 1) { | |
359 | pathInSi = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit)); | |
360 | PmodRecP = itsHit->GetParticle()->P(); | |
361 | PmodRecP *= 1.0e+3; | |
362 | // cout<<"path in Si="<<pathInSi<<endl; | |
363 | ||
364 | Float_t Signal = signalRP*(300/pathInSi)/38; | |
365 | if(Signal<0.5) { | |
366 | //cout<<"track,partcode,pmod,hitprim,hitstat,Signal,dEn ="<<trackHit<<","<<partcode<<","<<pmod<<","<<hitprim<<","<<hitstat<<","<<Signal<<","<<dEn<<endl; | |
367 | } | |
368 | } | |
369 | } // hit loop | |
370 | ||
371 | if(RecPointPrim == 1) { | |
372 | ||
373 | //cout<<" SDD/SSD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl; | |
374 | ||
375 | Int_t xpartcode=gAlice->Particle(TrackRecPoint)->GetPdgCode(); | |
376 | ||
377 | if(mod < first2) { // SDD | |
378 | signalRP *= (280./pathInSi); // 280 microns of Si thickness | |
379 | signalSDD->Fill(signalRP); | |
380 | signalSDDmip->Fill(signalRP/280.); // ADC/280 -> mips | |
381 | InPID(pid,(Int_t)nev,(Int_t)TrackRecPoint, signalRP/280. , | |
382 | (Float_t)PmodRecP,(Int_t)xpartcode); | |
383 | //cout<<" SDD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl; | |
384 | if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathInSi in SDD ="<<pathInSi<<endl; | |
385 | ||
386 | //if(signalRP/280 < 5) cout<<" small signalSDDmip, path ="<<signalRP/280<<","<<pathInSDD<<endl; | |
387 | }else{ // SSD | |
388 | signalRP *= (300/pathInSi); // 300 microns of Si thickness | |
389 | //cout<<" SSD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl; | |
390 | signalSSD->Fill(signalRP); | |
391 | signalSSDmip->Fill(signalRP/38.); // ADC/38 -> mips | |
392 | InPID(pid,(Int_t)nev,(Int_t)TrackRecPoint, signalRP/38. , | |
393 | (Float_t)PmodRecP,(Int_t)xpartcode); | |
394 | if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathInSi in SSD ="<<pathInSi<<endl; | |
395 | ||
396 | if(signalRP/38 < 0.5) cout<<" small signalSSD (mip), path, module ="<<signalRP/38<<","<<pathInSi<<mod<<endl; | |
397 | } | |
398 | } // primery particle | |
399 | } // pnt, recpoint loop | |
400 | ||
401 | // dEdX hit analysis for the SDD/SSD | |
402 | ||
403 | Int_t track = -100; | |
404 | Int_t trackPrev = -100; | |
405 | parent = -100; | |
406 | Int_t parentPrev = -100; | |
407 | Int_t flagprim = 0; | |
408 | Int_t TrackPart; | |
409 | epart = 0; | |
410 | pathInSi = 1.0e+7; | |
411 | Float_t PmodPart; | |
412 | zhit0 = 1.0e+7; | |
413 | xhit0 = 1.0e+7; | |
414 | yhit0 = 1.0e+7; | |
415 | ||
416 | /* | |
417 | if(mod <= 259) { | |
418 | cout<<" SDD hits: nhits ="<<nhits<<endl; | |
419 | }else{ | |
420 | cout<<" SSD hits: nhits ="<<nhits<<endl; | |
421 | } | |
422 | */ | |
423 | ||
424 | for (hit=0;hit<nhits;hit++) { | |
425 | ||
426 | itsHit = (AliITShit*)Mod->GetHit(hit); | |
427 | if(hit>0) itsHitPrev = (AliITShit*)Mod->GetHit(hit-1); | |
428 | hitprim = 0; | |
429 | Int_t hitprimPrev = 0; | |
430 | track = itsHit->GetTrack(); | |
431 | if(hit > 0) Int_t trackPrev = itsHitPrev->GetTrack(); | |
432 | dray = 0; | |
433 | hitstat = itsHit->GetTrackStatus(); | |
434 | ||
435 | zhit = 10000*itsHit->GetZL(); | |
436 | xhit = 10000*itsHit->GetXL(); | |
437 | yhit = 10000*itsHit->GetYL(); | |
438 | Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV | |
439 | ||
440 | // cout<<"New hit,lay,hitstat,yhit,ehit ="<<hit<<","<<layer<<","<<hitstat<<","<<yhit<<","<<ehit<<endl; | |
441 | // cout<<"befor: track,trackPrev ="<<track<<","<<trackPrev<<endl; | |
442 | ||
443 | parent = itsHit->GetParticle()->GetFirstMother(); | |
444 | if(hit>0) Int_t parentPrev = itsHitPrev->GetParticle()->GetFirstMother(); | |
445 | partcode = itsHit->GetParticle()->GetPdgCode(); | |
446 | ||
447 | // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+ | |
448 | // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda | |
449 | ||
450 | px = itsHit->GetPXL(); | |
451 | py = itsHit->GetPYL(); | |
452 | pz = itsHit->GetPZL(); | |
453 | pmod = 1000*sqrt(px*px+py*py+pz*pz); | |
454 | ||
455 | // cout<<"partcode,pmod,parent,parentPrev,hitprim,flagprim ="<<partcode<<","<<pmod<<","<<parent<<","<<parentPrev<<","<<hitprim<<","<<flagprim<<endl; | |
456 | ||
457 | ||
458 | if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- | |
459 | // at p < 6 MeV/c | |
460 | ||
461 | ||
462 | if(parent < 0 && parent > -90) hitprim += 1; | |
463 | if(parentPrev < 0 && parentPrev >-90) hitprimPrev += 1; | |
464 | // hitprim=1 for the primery particles | |
465 | // if(hit==0&&hitprim==1&&hitstat==66) { | |
466 | if(hitprim==1&&hitstat==66) { | |
467 | zhit0 = zhit; | |
468 | xhit0 = xhit; | |
469 | yhit0 = yhit; | |
470 | } | |
471 | ||
472 | if((hitstat == 68 || hitstat == 33) && hitprim == 1) { | |
473 | pathInSi = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit)); | |
474 | TrackPart = track; | |
475 | PmodPart = itsHit->GetParticle()->P(); | |
476 | PmodPart *= 1.0e+3; | |
477 | // cout<<"xhit0,yhit0,zhit0,pathInSi ="<<xhit0<<","<<yhit0<<","<<zhit0<<","<<pathInSi<<endl; | |
478 | } | |
479 | ||
480 | if(hitprim == 0) track = parent; | |
481 | if(hitprimPrev == 0) trackPrev = parentPrev; | |
482 | // cout<<"after: track,trackPrev ="<<track<<","<<trackPrev<<endl; | |
483 | ||
484 | if(hit > 0 && track == trackPrev) epart += ehit; | |
485 | if((hit > 0 && track == trackPrev) && (hitprim==1)) flagprim +=1; | |
486 | ||
487 | // cout<<"new hitprim, flagprim ="<<hitprim<<","<<flagprim<<endl; | |
488 | ||
489 | if((hit > 0 && track != trackPrev) || (hit == nhits-1)) { | |
490 | if(flagprim > 0) { | |
491 | if(layer > 2 && layer < 5) { | |
492 | if(epart < 40) cout<<" small SDD dedx ="<<epart<<endl; | |
493 | if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathSDD ="<<pathInSi<<endl; | |
494 | if(epart > 40 && pathInSi > 100 && pathInSi < 1.0e+5) { | |
495 | dedxSDD->Fill(epart); | |
496 | epart *= (280/pathInSi); | |
497 | dedxSDDnorm->Fill(epart); | |
498 | dedxSDDmip->Fill(epart/75); | |
499 | //InPID(pid,(Int_t)nev,(Int_t)track, epart/75. ,(Float_t)pmod,(Int_t)partcode); | |
500 | pathSDD->Fill(pathInSi); | |
501 | // cout<<"SDD hit: lay,lad,det,track,pmod,epart,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackPart<<","<<PmodPart<<","<<epart<<","<<pathInSi<<endl; | |
502 | } | |
503 | } | |
504 | if(layer > 4) { | |
505 | if(epart < 40) cout<<" small SSD dedx ="<<epart<<endl; | |
506 | if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathSSD ="<<pathInSi<<endl; | |
507 | if(epart > 40 && pathInSi > 100 && pathInSi < 1.0e+5) { | |
508 | dedxSSD->Fill(epart); | |
509 | epart *= (280/pathInSi); | |
510 | dedxSSDnorm->Fill(epart); | |
511 | dedxSSDmip->Fill(epart/80); | |
512 | //InPID(pid,(Int_t)nev,(Int_t)track, epart/80. ,(Float_t)pmod,(Int_t)partcode); | |
513 | pathSSD->Fill(pathInSi); | |
514 | //cout<<"SSD hit: lay,lad,det,track,pmod,epart,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackPart<<","<<PmodPart<<","<<epart<<","<<pathInSi<<endl; | |
515 | } | |
516 | } | |
517 | flagprim=0; | |
518 | epart = 0; | |
519 | }else{ | |
520 | //cout<<"No!, epart for the secondary"<<endl; | |
521 | epart = 0; | |
522 | } | |
523 | } | |
524 | } // SDD/SSD hit loop | |
525 | } // module loop | |
526 | //cout<<"negtrSDD, negtrSSD ="<<negtrSDD<<","<<negtrSSD<<endl; | |
527 | } // event loop | |
528 | ||
529 | cout<<endl; | |
530 | ||
531 | TFile fhistos("dEdX_his.root","RECREATE"); | |
532 | ||
533 | dedxSDD->Write(); | |
534 | dedxSDDnorm->Write(); | |
535 | dedxSDDmip->Write(); | |
536 | dedxSSD->Write(); | |
537 | dedxSSDnorm->Write(); | |
538 | dedxSSDmip->Write(); | |
539 | signalSDD->Write(); | |
540 | signalSDDmip->Write(); | |
541 | signalSSD->Write(); | |
542 | signalSSDmip->Write(); | |
543 | pathSDD->Write(); | |
544 | pathSSD->Write(); | |
545 | ||
546 | fhistos.Close(); | |
547 | cout<<"!!! Histogramms and ntuples were written"<<endl; | |
548 | // ------------------------------------------------------------ | |
549 | ||
550 | TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700); | |
551 | c1->Divide(2,2); | |
552 | ||
553 | ||
554 | // c1->cd(1); | |
555 | // gPad->SetFillColor(33); | |
556 | // signalSDD->SetFillColor(42); | |
557 | // signalSDD->Draw(); | |
558 | ||
559 | // c1->cd(2); | |
560 | // gPad->SetFillColor(33); | |
561 | // signalSDDmip->SetFillColor(42); | |
562 | // signalSDDmip->Draw(); | |
563 | ||
564 | // c1->cd(3); | |
565 | // gPad->SetFillColor(33); | |
566 | // dedxSDDnorm->SetFillColor(46); | |
567 | // dedxSDDnorm->Draw(); | |
568 | ||
569 | // c1->cd(4); | |
570 | // gPad->SetFillColor(33); | |
571 | // dedxSDDmip->SetFillColor(46); | |
572 | // dedxSDDmip->Draw(); | |
573 | ||
574 | c1->cd(1); | |
575 | gPad->SetFillColor(33); | |
576 | signalSSD->SetFillColor(42); | |
577 | signalSSD->Draw(); | |
578 | ||
579 | c1->cd(2); | |
580 | gPad->SetFillColor(33); | |
581 | signalSSDmip->SetFillColor(42); | |
582 | signalSSDmip->Draw(); | |
583 | ||
584 | c1->cd(3); | |
585 | gPad->SetFillColor(33); | |
586 | dedxSSDnorm->SetFillColor(46); | |
587 | dedxSSDnorm->Draw(); | |
588 | ||
589 | c1->cd(4); | |
590 | gPad->SetFillColor(33); | |
591 | dedxSSDmip->SetFillColor(46); | |
592 | dedxSSDmip->Draw(); | |
593 | ||
594 | ||
595 | /* | |
596 | c1->cd(1); | |
597 | gPad->SetFillColor(33); | |
598 | dedxSDD->SetFillColor(42); | |
599 | dedxSDD->Draw(); | |
600 | ||
601 | c1->cd(2); | |
602 | gPad->SetFillColor(33); | |
603 | pathSDD->SetFillColor(42); | |
604 | pathSDD->Draw(); | |
605 | ||
606 | c1->cd(3); | |
607 | gPad->SetFillColor(33); | |
608 | dedxSSD->SetFillColor(46); | |
609 | dedxSSD->Draw(); | |
610 | ||
611 | c1->cd(4); | |
612 | gPad->SetFillColor(33); | |
613 | pathSSD->SetFillColor(46); | |
614 | pathSSD->Draw(); | |
615 | */ | |
616 | ||
617 | // if(!pid)pid->Tab(); | |
618 | cout<<"END test for clusters and hits "<<endl; | |
619 | ||
620 | } | |
621 | //============================ InPID() ========================================= | |
622 | int totpid; | |
623 | void InPID (AliITSPid *pid=0,int nev, | |
624 | Int_t track, | |
625 | Float_t signal, | |
626 | Float_t pmom, | |
627 | Int_t pcode) | |
628 | { | |
629 | // Includes recp in PID table. | |
630 | if(pid){ | |
631 | //cout<<" In pid: track,sig,pmod,pcode="<<track<<","<<signal<<","<<pmom<<","<<pcode<<endl; | |
632 | if( abs(pcode)==321 || abs(pcode)==211 ||abs(pcode)==2212 ||abs(pcode)==11 ) | |
633 | { | |
634 | pid->SetEdep(10000*nev+track,signal); | |
635 | pid->SetPmom(10000*nev+track,pmom); | |
636 | pid->SetPcod(10000*nev+track,abs(pcode)); | |
637 | } | |
638 | } | |
639 | totpid++; | |
640 | } | |
641 | //======================= | |
642 | ||
643 | ||
644 | ||
645 | ||
646 | ||
647 |