Partical update of SetDefault.
[u/mrichter/AliRoot.git] / ITS / SSDrecpointTest.C
CommitLineData
6b8f55ce 1void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
2 //void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=999)
3{
4/////////////////////////////////////////////////////////////////////////
5// This macro is a small example of a ROOT macro
6// illustrating how to read the output of GALICE
7// and fill some histograms.
8//
9// Root > .L anal.C //this loads the macro in memory
10// Root > anal(); //by default process first event
11// Root > anal(2); //process third event
12//Begin_Html
13/*
14<img src="gif/anal.gif">
15*/
16//End_Html
17/////////////////////////////////////////////////////////////////////////
18
19// Dynamically link some shared libs
20
21 if (gClassTable->GetID("AliRun") < 0) {
22 gROOT->LoadMacro("loadlibs.C");
23 loadlibs();
24 }
25
26// Connect the Root Galice file containing Geometry, Kine and Hits
27 TString *str = new TString("galice.root");
28 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data());
29 if (!file) file = new TFile(str->Data(),"UPDATE");
30
31// Get AliRun object from file or create it if not on file
32 // if (!gAlice) {
33 gAlice = (AliRun*)file->Get("gAlice");
34 if (gAlice) printf("AliRun object found on file\n");
35 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
36 //}
37
38
39 // -------------- Create ntuples --------------------
40
41 // ntuple structures:
42
43
44 struct {
45 Int_t lay;
46 Int_t nxP;
47 Int_t nxN;
48 Int_t hitprim;
49 Int_t partcode;
50 Float_t x;
51 Float_t z;
52 Float_t dx;
53 Float_t dz;
54 Float_t pmod;
55 } ntuple_st;
56
57 struct {
58 Int_t lay;
59 Int_t lad;
60 Int_t det;
61 Int_t nxP;
62 Int_t nxN;
63 Int_t noverlaps;
64 Int_t noverprim;
65 Float_t qclP;
66 Float_t qclN;
67 Float_t qrec;
68 Float_t dx;
69 Float_t dz;
70 } ntuple1_st;
71
72 struct {
73 Int_t nxP;
74 Int_t nxN;
75 Float_t x;
76 Float_t z;
77 } ntuple2_st;
78
79 ntuple = new TTree("ntuple","Demo ntuple");
80 ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
81 ntuple->Branch("nxP",&ntuple_st.nxP,"nxP/I");
82 ntuple->Branch("nxN",&ntuple_st.nxN,"nxN/I");
83 ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
84 ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
85 ntuple->Branch("x",&ntuple_st.x,"x/F");
86 ntuple->Branch("z",&ntuple_st.z,"z/F");
87 ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
88 ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
89 ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
90
91 ntuple1 = new TTree("ntuple1","Demo ntuple1");
92 ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
93 ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
94 ntuple1->Branch("det",&ntuple1_st.det,"det/I");
95 ntuple1->Branch("nxP",&ntuple1_st.nxP,"nxP/I");
96 ntuple1->Branch("nxN",&ntuple1_st.nxN,"nxN/I");
97 ntuple1->Branch("qclP",&ntuple1_st.qclP,"qclP/F");
98 ntuple1->Branch("qclN",&ntuple1_st.qclN,"qclN/F");
99 ntuple1->Branch("qrec",&ntuple1_st.qrec,"qrec/F");
100 ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
101 ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
102 ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
103 ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
104
105 ntuple2 = new TTree("ntuple2","Demo ntuple2");
106 ntuple2->Branch("nxP",&ntuple2_st.nxP,"nxP/I");
107 ntuple2->Branch("nxN",&ntuple2_st.nxN,"nxN/I");
108 ntuple2->Branch("x",&ntuple2_st.x,"x/F");
109 ntuple2->Branch("z",&ntuple2_st.z,"z/F");
110
111
112 // Create Histogramms
113
114 TH1F *NxP5 = new TH1F("NxP5","P cluster size for layer 5",20,0.,20.);
115 TH1F *NxN5 = new TH1F("NxN5","N cluster size for layer 5",20,0.,20.);
116 TH1F *NxP6 = new TH1F("NxP6","P cluster size for layer 6",20,0.,20.);
117 TH1F *NxN6 = new TH1F("NxN6","N cluster size for layer 6",20,0.,20.);
118
119 TH1F *Xres5 = new TH1F("Xres5","Xrec and Xgen difference (micr) for layers 5",100,-200.,200.);
120 TH1F *Xres6 = new TH1F("Xres6","Xrec and Xgen difference (micr) for layers 6",100,-200.,200.);
121 TH1F *Zres5 = new TH1F("Zres5","Zrec and Zgen difference (micr) for layers 5",100,-8000.,8000.);
122 TH1F *Zres6 = new TH1F("Zres6","Zrec and Zgen difference (micr) for layers 6",100,-8000.,8000.);
123 TH1F *Path5 = new TH1F("Path5","Path length in Si",100,0.,600.);
124 TH1F *Path6 = new TH1F("Path6","Path length in Si",100,0.,600.);
125 TH1F *dEdX = new TH1F("dEdX","dEdX (KeV)",100,0.,500.);
126 TH2F *adcPadcN5all = new TH2F("adcPadcN5all","adcP/N correlation for lay5",100,0.,200.,100,0.,200.);
127 TH2F *adcPadcN6all = new TH2F("adcPadcN6all","adcP/N correlation for lay6",100,0.,200.,100,0.,200.);
128 TH2F *adcPadcN5cut = new TH2F("adcPadcN5cut","adcP/N correlation for lay5 and cut of P-N signas",100,0.,200.,100,0.,200.);
129 TH2F *adcPadcN6cut = new TH2F("adcPadcN6cut","adcP/N correlation for lay6 and cut of P-N signals",100,0.,200.,100,0.,200.);
130
131
132 AliITS *ITS = (AliITS*) gAlice->GetModule("ITS");
133 if (!ITS) { cout << "no ITS" << endl; return; }
134
3f0d03a1 135 //AliITSgeom *aliitsgeo = ITS->GetITSgeom();
136 AliITSgeom *geom = ITS->GetITSgeom();
6b8f55ce 137
138 //Int_t cp[8]={0,0,0,0,0,0,0,0};
139
140 cout << "SSD" << endl;
141
142 AliITSDetType *iDetType=ITS->DetType(2);
143 AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
144 AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
145 //res2->SetSigmaSpread(3.,2.);
146 AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
147 ITS->SetSimulationModel(2,sim2);
148
149 TClonesArray *dig2 = ITS->DigitsAddress(2);
150 TClonesArray *recp2 = ITS->ClustersAddress(2);
151 // AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2,recp2);
152 AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
153 ITS->SetReconstructionModel(2,rec2);
154 // test
155 printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz());
156 printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
157
158
159//
160// Loop over events
161//
162
163
164 Int_t Nh=0;
165 Int_t Nh1=0;
166 for (int nev=0; nev<= evNumber2; nev++) {
167 Int_t nparticles = 0;
168 nparticles = gAlice->GetEvent(nev);
169 cout << "nev " << nev <<endl;
170 cout << "nparticles " << nparticles <<endl;
171 if (nev < evNumber1) continue;
172 if (nparticles <= 0) return;
173
174 AliITShit *itsHit;
175 AliITSRecPoint *itsPnt = 0;
176 AliITSRawClusterSSD *itsClu = 0;
177
178 // Get Hit, Cluster & Recpoints Tree Pointers
179
180 TTree *TH = gAlice->TreeH();
181 Int_t nenthit=TH->GetEntries();
182 printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
183
184 ITS->GetTreeC(nev);
185 TTree *TC=ITS->TreeC();
186 Int_t nentclu=TC->GetEntries();
187 printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
188
189 TTree *TR = gAlice->TreeR();
190 Int_t nentrec=TR->GetEntries();
191 printf("Found %d entries in the RecPoints tree\n",nentrec);
192
193 // Get Pointers to Clusters & Recpoints TClonesArrays
194
195 TClonesArray *ITSclu = ITS->ClustersAddress(2);
196 printf ("ITSclu %p \n",ITSclu);
197 TClonesArray *ITSrec = ITS->RecPoints();
198 printf ("ITSrec %p \n",ITSrec);
199
200 // check recpoints
201
202 //Int_t nbytes = 0;
203 Int_t totpoints = 0;
204 Int_t totclust = 0;
205
206 // check hits
207
208 Int_t nmodules=0;
3f0d03a1 209 Int_t mod;
6b8f55ce 210
211 ITS->InitModules(-1,nmodules);
212 ITS->FillModules(nev,0,nmodules,"","");
213
214 TObjArray *fITSmodules = ITS->GetModules();
215
3f0d03a1 216 Int_t first0 = geom->GetStartDet(0); // SPD
217 Int_t last0 = geom->GetLastDet(0); // SPD
218 Int_t first1 = geom->GetStartDet(1); // SDD
219 Int_t last1 = geom->GetLastDet(1); // SDD
220 Int_t first2 = geom->GetStartDet(2); // SSD
221 Int_t last2 = geom->GetLastDet(2); // SSD
222
223 // For the SPD: first0 = 0, last0 = 239 (240 modules);
224 // for the SDD: first1 = 240, last1 = 499 (260 modules);
225 // for the SSD: first2 = 500, last2 = 2269 (1770 modules).
226
227 printf("det type %d first0, last0 %d %d \n",0,first0,last0);
228 printf("det type %d first1, last1 %d %d \n",1,first1,last1);
229 printf("det type %d first2, last2 %d %d \n",2,first2,last2);
6b8f55ce 230
3f0d03a1 231 // module loop for the SSD
232 for (mod=first2; mod<last2+1; mod++) {
6b8f55ce 233
3f0d03a1 234 TTree *TR = gAlice->TreeR();
6b8f55ce 235 Int_t nentrec=TR->GetEntries();
236 //printf("Found %d entries in the RecPoints tree\n",nentrec);
237
238 //cout << "CLUSTERS: reset" << endl;
239 ITS->ResetClusters();
240 //cout << "CLUSTERS: get" << endl;
241 TC->GetEvent(mod);
242 //cout << "RECPOINTS: reset" << endl;
243 ITS->ResetRecPoints();
244 //cout << "RECPOINTS: get" << endl;
3f0d03a1 245 //TR->GetEvent(mod+1); // for the V3.04 AliRoot
246 TR->GetEvent(mod); // for the V3.05 AliRoot
6b8f55ce 247
248 Int_t nrecp = ITSrec->GetEntries();
249 totpoints += nrecp;
250 if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
251 if (!nrecp) continue;
252 Int_t nclusters = ITSclu->GetEntries();
253 totclust += nclusters;
6b8f55ce 254 //if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
255
256 AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
257 // printf("Mod: %X\n",Mod);
258 Int_t nhits = Mod->GetNhits();
259 Float_t epart = 0;
260 cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
261
262 // ---------------- cluster/hit analysis ---------------------
263
264
265 Float_t pathInSSD = 300.;
266
267 // ---- Recpoint loop
268 for (Int_t pnt=0;pnt<nrecp;pnt++) {
269 itsPnt = (AliITSRecPoint*)ITSrec->At(pnt);
270 if(!itsPnt) continue;
271 itsClu = (AliITSRawClusterSSD*)ITSclu->At(pnt);
272 if(!itsClu) continue;
273
274 Int_t nxP = itsClu->fMultiplicity;
275 Int_t nxN = itsClu->fMultiplicityN;
276 Float_t qclP = itsClu->fSignalP; // in ADC
277 Float_t qclN = itsClu->fSignalN; // in ADC
278 //Float_t dq = qclP - qclN;
279 Float_t qcut = itsClu->fQErr; // abs(dq)/signal,
280 // where signal is
281 // max of qclP,qclN
282 Float_t xrec = 10000*itsPnt->GetX();
283 Float_t zrec = 10000*itsPnt->GetZ();
284 Float_t qrec = itsPnt->GetQ(); // in ADC, maximum from fSignalP/N
285 //Float_t dedx = itsPnt->GetdEdX(); // in KeV (ADC * 2.16)
286 Float_t dedx = itsPnt->fdEdX; // in KeV (ADC * 2.16)
287 Int_t ii = 0;
288 Int_t tr1 = itsPnt->GetLabel(ii);
289 Int_t ii = 1;
290 Int_t tr2 = itsPnt->GetLabel(ii);
291 Int_t ii = 2;
292 Int_t tr3 = itsPnt->GetLabel(ii);
293
294 // fill ntuple2
295 ntuple2_st.nxP = nxP;
296 ntuple2_st.nxN = nxN;
297 ntuple2_st.x = xrec/1000;
298 ntuple2_st.z = zrec/1000;
299
300 if(qcut < 0.18) ntuple2->Fill();
301
302
303 Int_t noverlaps = 0;
304 Int_t noverprim = 0;
305 Int_t flaghit = 0;
306 Float_t xhit0 = 1e+7;
307 Float_t yhit0 = 1e+7;
308 Float_t zhit0 = 1e+7;
309
310 // Hit loop
311 for (Int_t hit=0;hit<nhits;hit++) {
312
313 itsHit = (AliITShit*)Mod->GetHit(hit);
314
315 Int_t flagtrack = 0;
316 Int_t hitlayer = itsHit->GetLayer();
317 Int_t hitladder= itsHit->GetLadder();
318 Int_t hitdet= itsHit->GetDetector();
319
320 Int_t track = itsHit->GetTrack();
321 Int_t dray = 0;
322 Int_t hitstat = itsHit->GetTrackStatus();
323
324 Float_t zhit = 10000*itsHit->GetZL();
325 Float_t xhit = 10000*itsHit->GetXL();
326 Float_t yhit = 10000*itsHit->GetYL();
327 Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
328
329 Int_t parent = itsHit->GetParticle()->GetFirstMother();
330 Int_t partcode = itsHit->GetParticle()->GetPdgCode();
331
332 // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
333 // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
334
335 Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
336 // vertex
337 pmod *= 1.0e+3;
338
339 if(hitstat == 66 && yhit < -146.) { // entering hit
340 xhit0 = xhit;
341 yhit0 = yhit;
342 zhit0 = zhit;
343 }
344
345 if(hitstat == 66) continue; // Take the not entering hits only
346
347 if(xhit0 > 9e+6 || zhit0 > 9e+6 || yhit0 > 9e+6) {
348 //cout<<"default xhit0,zhit0,yhit0 ="<<xhit0<<","<<zhit0<<","<<yhit0<<endl;
349 continue;
350 }
351
352
353
354 // Consider the hits only with the track number equaled to one
355 // of the recpoint
356 if(track == tr1) flagtrack = 1;
357
358 if(flagtrack == 1) { // the hit corresponds to the recpoint
359
360 flaghit = 1;
361
362 //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
363 //Float_t py = itsHit->GetPYL();
364 //Float_t pz = itsHit->GetPZL();
365
366 Int_t hitprim = 0;
367
368 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
369 // at p < 6 MeV/c
370
371 if((hitstat == 68 || hitstat == 33) && dray == 0) noverlaps=noverlaps + 1;
372 // overlapps for all hits but
373 // not for delta ray which
374 // also went out from the
375 // detector and returned
376 // again
377
378
379 // x,z resolution colculation
380 if((hitstat == 68 || hitsat == 33) && parent < 0) {
381 Float_t xmed = (xhit + xhit0)/2;
382 Float_t zmed = (zhit + zhit0)/2;
383 Float_t xdif = xmed - xrec;
384 Float_t zdif = zmed - zrec;
385
386 hitprim = 1; // hitprim=1 for the primery particles
387
388 noverprim += 1;
389
390 pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
391
392
393 // fill ntuple
394 ntuple_st.lay = hitlayer;
395 ntuple_st.nxP = nxP;
396 ntuple_st.nxN = nxN;
397 ntuple_st.hitprim = hitprim;
398 ntuple_st.partcode = partcode;
399 ntuple_st.x = xrec/1000;
400 ntuple_st.z = zrec/1000;
401 ntuple_st.dx = xdif;
402 ntuple_st.dz = zdif;
403 ntuple_st.pmod = pmod;
404
405 if(qcut < 0.18) ntuple->Fill();
406
407 if(hitlayer == 5 && qcut < 0.18) {
408 Xres5->Fill(xdif);
409 Zres5->Fill(zdif);
410 Path5->Fill(pathInSSD);
411 }
412 if(hitlayer == 6 && qcut < 0.18) {
413 Xres6->Fill(xdif);
414 Zres6->Fill(zdif);
415 Path6->Fill(pathInSSD);
416 }
417 } // hitstat 68/33
418 } else { // non correspondent hit
419 xhit0 = 1e+7;
420 zhit0 = 1e+7;
421 } // end of hit-recpoint correspondence
422 } // hit loop
423
424 if(flaghit == 1) {
6b8f55ce 425
426 if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
427 // delta rays only
428
429 // fill ntuple1
430 ntuple1_st.lay = hitlayer;
431 ntuple1_st.lad = hitladder;
432 ntuple1_st.det = hitdet;
433 ntuple1_st.nxP = nxP;
434 ntuple1_st.nxN = nxN;
435 ntuple1_st.qclP = qclP*300/pathInSSD;
436 ntuple1_st.qclN = qclN*300/pathInSSD;
437 ntuple1_st.qrec = qrec*300/pathInSSD;
438 ntuple1_st.dx = xdif;
439 ntuple1_st.dz = zdif;
440 noverlaps -= 1;
441 noverprim -= 1;
442 ntuple1_st.noverlaps = noverlaps;
443 ntuple1_st.noverprim = noverprim;
444
445 if(qcut < 0.18) ntuple1->Fill();
446
447 Float_t de = dedx*300./pathInSSD;
448 dEdX->Fill(de);
449 if(hitlayer == 5 ) {
450 adcPadcN5all->Fill(qclP,qclN);
451 }
452 if(hitlayer == 6 ) {
453 adcPadcN6all->Fill(qclP,qclN);
454 }
455 if(hitlayer == 5 && qcut < 0.18) {
456 adcPadcN5cut->Fill(qclP,qclN);
457 NxP5->Fill(nxP);
458 NxN5->Fill(nxN);
459 }
460 if(hitlayer == 6 && qcut < 0.18) {
461 adcPadcN6cut->Fill(qclP,qclN);
462 NxP6->Fill(nxP);
463 NxN6->Fill(nxN);
464 }
465 } // flaghit = 1
466 } //b.b. recpoint loop
467 } //b.b. module loop
468 } //b.b. evnt loop
469
470 TFile fhistos("SSD_his.root","RECREATE");
471
472 ntuple->Write();
473 ntuple1->Write();
474 ntuple2->Write();
475 NxP5->Write();
476 NxN5->Write();
477 NxP6->Write();
478 NxN6->Write();
479 Xres5->Write();
480 Zres5->Write();
481 Xres6->Write();
482 Zres6->Write();
483 Path5->Write();
484 Path6->Write();
485 adcPadcN5all->Write();
486 adcPadcN6all->Write();
487 adcPadcN5cut->Write();
488 adcPadcN6cut->Write();
489 dEdX->Write();
490
491 fhistos.Close();
492
493 cout<<"!!! Histogramms and ntuples were written"<<endl;
494
495 TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
496 c1->Divide(2,2);
497 c1->cd(1);
498 gPad->SetFillColor(33);
499 Xres5->SetFillColor(42);
500 Xres5->Draw();
501 c1->cd(2);
502 gPad->SetFillColor(33);
503 Zres5->SetFillColor(46);
504 Zres5->Draw();
505 c1->cd(3);
506 gPad->SetFillColor(33);
507 Xres6->SetFillColor(42);
508 Xres6->Draw();
509 c1->cd(4);
510 gPad->SetFillColor(33);
511 Zres6->SetFillColor(46);
512 Zres6->Draw();
513
514 cout<<"END test for clusters and hits "<<endl;
515
516}
517
518
519
520
521
522