Test macros for SSD reconstruction
[u/mrichter/AliRoot.git] / ITS / SSDrecpointTest.C
1 void 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    
135    AliITSgeom *aliitsgeo = ITS->GetITSgeom();
136
137    //Int_t cp[8]={0,0,0,0,0,0,0,0};
138
139    cout << "SSD" << endl;
140
141    AliITSDetType *iDetType=ITS->DetType(2);
142    AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
143    AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
144    //res2->SetSigmaSpread(3.,2.);
145    AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
146    ITS->SetSimulationModel(2,sim2);
147
148    TClonesArray *dig2  = ITS->DigitsAddress(2);
149    TClonesArray *recp2  = ITS->ClustersAddress(2);
150    //   AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2,recp2);
151    AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
152    ITS->SetReconstructionModel(2,rec2);
153    // test
154    printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz());
155    printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
156
157    
158 //
159 //   Loop over events
160 //
161
162
163    Int_t Nh=0;
164    Int_t Nh1=0;
165    for (int nev=0; nev<= evNumber2; nev++) {
166      Int_t nparticles = 0;
167      nparticles = gAlice->GetEvent(nev);
168      cout << "nev         " << nev <<endl;
169      cout << "nparticles  " << nparticles <<endl;
170      if (nev < evNumber1) continue;
171      if (nparticles <= 0) return;
172      
173      AliITShit *itsHit;
174      AliITSRecPoint *itsPnt = 0;
175      AliITSRawClusterSSD *itsClu = 0;
176      
177      // Get Hit, Cluster & Recpoints Tree Pointers
178
179      TTree *TH = gAlice->TreeH();
180      Int_t nenthit=TH->GetEntries();
181      printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
182
183      ITS->GetTreeC(nev);
184      TTree *TC=ITS->TreeC();
185      Int_t nentclu=TC->GetEntries();
186      printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
187
188      TTree *TR = gAlice->TreeR();
189      Int_t nentrec=TR->GetEntries();
190      printf("Found %d entries in the RecPoints tree\n",nentrec);
191
192      // Get Pointers to Clusters & Recpoints TClonesArrays
193
194      TClonesArray *ITSclu  = ITS->ClustersAddress(2); 
195      printf ("ITSclu %p \n",ITSclu);
196      TClonesArray *ITSrec  = ITS->RecPoints(); 
197      printf ("ITSrec %p \n",ITSrec);
198
199      // check recpoints
200
201      //Int_t nbytes = 0;
202      Int_t totpoints = 0;
203      Int_t totclust = 0;
204
205      // check hits
206      
207      Int_t nmodules=0;
208      
209      ITS->InitModules(-1,nmodules); 
210      ITS->FillModules(nev,0,nmodules,"","");
211      
212      TObjArray *fITSmodules = ITS->GetModules();
213      
214      Int_t first1 = aliitsgeo->GetStartDet(2); // SSD
215      Int_t last1 = aliitsgeo->GetLastDet(2);   // SSD
216      Int_t first = aliitsgeo->GetStartDet(1);  // SDD
217      Int_t last = aliitsgeo->GetLastDet(1);    // SDD
218      Int_t first0 = aliitsgeo->GetStartDet(0);  // SPD
219      Int_t last0 = aliitsgeo->GetLastDet(0);    // SPD
220
221      //  For the SPD: first0 = 0, last0 = 239, 240 modules;  
222      //  for the SDD: first1 = 240, last1 = 499, 260 modules;  
223      //  for the SPD: first2 = 500, last2 = 2269, 1770 modules;
224        
225      printf("det type %d first, last %d %d \n",0,first0,last0);
226      printf("det type %d first, last %d %d \n",1,first,last);
227      printf("det type %d first1, last1 %d %d \n",2,first1,last1);
228
229      // module loop
230      for (Int_t mod=first1; mod<last1; mod++) {
231
232               TTree *TR = gAlice->TreeR();
233        Int_t nentrec=TR->GetEntries();
234        //printf("Found %d entries in the RecPoints tree\n",nentrec);
235       
236               //cout << "CLUSTERS: reset" << endl;
237        ITS->ResetClusters();
238        //cout << "CLUSTERS: get" << endl;
239        TC->GetEvent(mod);
240        //cout << "RECPOINTS: reset" << endl;
241        ITS->ResetRecPoints();
242        //cout << "RECPOINTS: get" << endl;
243        //TR->GetEvent(mod+1);
244        TR->GetEvent(mod);
245
246        Int_t nrecp = ITSrec->GetEntries();
247        totpoints += nrecp;
248        if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
249        if (!nrecp) continue;
250        Int_t nclusters = ITSclu->GetEntries();
251        totclust += nclusters;
252        if (nclusters) cout<<"ncluster ="<<nclusters<<endl;
253        //if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
254        
255        AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
256        //       printf("Mod: %X\n",Mod);
257        Int_t nhits = Mod->GetNhits();
258        Float_t epart = 0;
259        cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
260
261        // ---------------- cluster/hit analysis ---------------------
262
263
264      Float_t pathInSSD = 300.;
265
266        // ---- Recpoint loop
267        for (Int_t pnt=0;pnt<nrecp;pnt++) {
268          itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
269          if(!itsPnt) continue;
270          itsClu  = (AliITSRawClusterSSD*)ITSclu->At(pnt);
271          if(!itsClu) continue;
272
273          Int_t nxP = itsClu->fMultiplicity;
274          Int_t nxN = itsClu->fMultiplicityN;
275          Float_t qclP = itsClu->fSignalP;     // in ADC
276          Float_t qclN = itsClu->fSignalN;     // in ADC
277          //Float_t dq = qclP - qclN;
278          Float_t qcut = itsClu->fQErr;        // abs(dq)/signal,
279                                               // where signal is
280                                               // max of qclP,qclN        
281          Float_t xrec = 10000*itsPnt->GetX();
282          Float_t zrec = 10000*itsPnt->GetZ();
283          Float_t qrec = itsPnt->GetQ();      // in ADC, maximum from fSignalP/N
284          //Float_t dedx = itsPnt->GetdEdX();   // in KeV (ADC * 2.16)
285          Float_t dedx = itsPnt->fdEdX;   // in KeV (ADC * 2.16)
286          Int_t ii = 0;
287          Int_t tr1 = itsPnt->GetLabel(ii);
288          Int_t ii = 1;
289          Int_t tr2 = itsPnt->GetLabel(ii);
290          Int_t ii = 2;
291          Int_t tr3 = itsPnt->GetLabel(ii);
292
293          // fill ntuple2
294              ntuple2_st.nxP = nxP;
295              ntuple2_st.nxN = nxN;
296              ntuple2_st.x = xrec/1000;
297              ntuple2_st.z = zrec/1000;
298
299              if(qcut < 0.18) ntuple2->Fill();
300
301
302           Int_t noverlaps = 0;
303           Int_t noverprim = 0;
304           Int_t flaghit = 0;
305           Float_t xhit0 = 1e+7;
306           Float_t yhit0 = 1e+7;
307           Float_t zhit0 = 1e+7;
308
309        // Hit loop
310         for (Int_t hit=0;hit<nhits;hit++) {
311
312          itsHit   = (AliITShit*)Mod->GetHit(hit);
313
314          Int_t flagtrack = 0;
315          Int_t hitlayer = itsHit->GetLayer();
316          Int_t hitladder= itsHit->GetLadder();
317          Int_t hitdet= itsHit->GetDetector();
318
319          Int_t track = itsHit->GetTrack();
320          Int_t dray = 0;
321          Int_t hitstat = itsHit->GetTrackStatus();
322
323           Float_t zhit = 10000*itsHit->GetZL();
324           Float_t xhit = 10000*itsHit->GetXL();
325           Float_t yhit = 10000*itsHit->GetYL();
326           Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV 
327
328            Int_t parent = itsHit->GetParticle()->GetFirstMother();
329            Int_t partcode = itsHit->GetParticle()->GetPdgCode();
330
331    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
332    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
333
334            Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
335                                                       // vertex
336            pmod *= 1.0e+3;
337
338           if(hitstat == 66 && yhit < -146.) {  // entering hit
339             xhit0 = xhit;
340             yhit0 = yhit;
341             zhit0 = zhit;
342           }
343
344           if(hitstat == 66) continue; // Take the not entering hits only 
345
346           if(xhit0 > 9e+6 || zhit0 > 9e+6 || yhit0 > 9e+6) {
347             //cout<<"default xhit0,zhit0,yhit0 ="<<xhit0<<","<<zhit0<<","<<yhit0<<endl;
348             continue;
349           }
350
351
352
353           // Consider the hits only with the track number equaled to one
354           // of the recpoint
355           if(track == tr1) flagtrack = 1;
356
357          if(flagtrack == 1) {     // the hit corresponds to the recpoint
358
359            flaghit = 1;
360
361            //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
362            //Float_t py = itsHit->GetPYL();
363            //Float_t pz = itsHit->GetPZL();
364
365          Int_t hitprim = 0;
366
367          if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
368                                                   // at p < 6 MeV/c
369
370          if((hitstat == 68 || hitstat == 33) && dray == 0)  noverlaps=noverlaps + 1;
371                                                   // overlapps for all hits but
372                                                   // not for delta ray which
373                                                   // also went out from the
374                                                   // detector and returned
375                                                   // again
376
377
378           // x,z resolution colculation
379           if((hitstat == 68 || hitsat == 33) && parent < 0)  {
380              Float_t xmed = (xhit + xhit0)/2;
381              Float_t zmed = (zhit + zhit0)/2;
382              Float_t xdif = xmed - xrec;
383              Float_t zdif = zmed - zrec;
384
385              hitprim = 1; // hitprim=1 for the primery particles
386
387              noverprim += 1;
388
389              pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
390
391
392          // fill ntuple
393              ntuple_st.lay = hitlayer;
394              ntuple_st.nxP = nxP;
395              ntuple_st.nxN = nxN;
396              ntuple_st.hitprim = hitprim;
397              ntuple_st.partcode = partcode;
398              ntuple_st.x = xrec/1000;
399              ntuple_st.z = zrec/1000;
400              ntuple_st.dx = xdif;
401              ntuple_st.dz = zdif;
402              ntuple_st.pmod = pmod;
403
404              if(qcut < 0.18) ntuple->Fill();
405
406             if(hitlayer == 5 && qcut < 0.18) {
407              Xres5->Fill(xdif);
408              Zres5->Fill(zdif);
409              Path5->Fill(pathInSSD);
410             }
411             if(hitlayer == 6 && qcut < 0.18) {
412              Xres6->Fill(xdif);
413              Zres6->Fill(zdif);
414              Path6->Fill(pathInSSD);
415             }
416           } // hitstat 68/33
417          } else {       // non correspondent hit
418           xhit0 = 1e+7;
419           zhit0 = 1e+7;
420          } // end of hit-recpoint correspondence
421         } // hit loop       
422
423         if(flaghit == 1) {
424           cout<<" flaghit Ok"<<endl;
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