]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/dEdXgeant.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated Linkdef and libTOF.pkg
[u/mrichter/AliRoot.git] / ITS / dEdXgeant.C
CommitLineData
23efe5f1 1
2Int_t max_modules=2269;
3
4void 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) {
270cout<<"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() =========================================
622int totpid;
623void 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