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