Release version of ITS code
[u/mrichter/AliRoot.git] / ITS / SPDclusterTest.C
1 #include "iostream.h"
2
3 void SPDclusterTest (Int_t evNumber1=0,Int_t evNumber2=0) 
4 {
5 /////////////////////////////////////////////////////////////////////////
6 //   This macro is a small example of a ROOT macro
7 //   illustrating how to read the output of GALICE
8 //   and do some analysis.
9 //   
10 /////////////////////////////////////////////////////////////////////////
11
12 // Dynamically link some shared libs
13
14    if (gClassTable->GetID("AliRun") < 0) {
15       gROOT->LoadMacro("loadlibs.C");
16       loadlibs();
17    }
18
19 // Connect the Root Galice file containing Geometry, Kine and Hits
20
21    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
22    if (!file) file = new TFile("galice.root");
23    file->ls();
24
25 // Get AliRun object from file or create it if not on file
26
27    if (!gAlice) {
28       gAlice = (AliRun*)file->Get("gAlice");
29       if (gAlice) printf("AliRun object found on file\n");
30       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
31    }
32  
33 //
34 //   Loop over events 
35 //
36    Int_t Nh=0;
37    Int_t Nh1=0;
38    for (int nev=0; nev<= evNumber2; nev++) {
39      Int_t nparticles = gAlice->GetEvent(nev);
40      cout << "nev         " << nev <<endl;
41      cout << "nparticles  " << nparticles <<endl;
42      if (nev < evNumber1) continue;
43      if (nparticles <= 0) return;
44
45      TTree *TH = gAlice->TreeH();
46      Int_t ntracks = TH->GetEntries();
47      cout<<"ntracks "<<ntracks<<endl;
48
49    Int_t nbytes = 0;
50
51    AliITSRawClusterSPD  *ITSclust;
52
53 // Get pointers to Alice detectors and Digits containers
54    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
55    TClonesArray *Particles = gAlice->Particles();
56
57    if (ITS) {
58      // fill modules with sorted by module hits
59      Int_t nmodules;
60      ITS->InitModules(-1,nmodules); 
61      ITS->FillModules(nev,-1,evNumber2,nmodules," "," ");
62      //get pointer to modules array
63      TObjArray *ITSmodules = ITS->GetModules();
64      AliITShit *itsHit;
65
66      // get the Tree for clusters
67      ITS->GetTreeC(nev);
68      TTree *TC=ITS->TreeC();
69      Int_t nent=TC->GetEntries();
70      printf("Found %d entries in the tree (must be one per module per event!)\n",nent);
71    
72      for (Int_t idettype=0;idettype<3;idettype++) {
73
74        TClonesArray *ITSclusters  = ITS->ClustersAddress(idettype);
75        //printf ("ITSclusters %p \n",ITSclusters);
76
77           if (idettype != 0) continue;
78
79
80           // ------------ Cluster and point analysis histogramms ------------
81
82           TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.);
83           TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.);
84           TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.);
85           TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.);
86           TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.);
87           TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.);
88           TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.);
89           TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.);
90
91           TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.);
92           TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.);
93           TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.);
94           TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.);
95
96
97           // -------------- Create ntuples --------------------
98
99           //  ntuple structures:
100
101           struct {
102             Int_t lay;
103             Int_t nx;
104             Int_t nz;
105             Int_t hitprim;
106             Int_t partcode;
107             Float_t dx;
108             Float_t dz;
109             Float_t pmod;
110           } ntuple_st;
111
112           struct {
113             Int_t lay;
114             Int_t lad;
115             Int_t det;
116             Int_t nx;
117             Int_t nz;
118             Int_t noverlaps;
119             Int_t noverprim;
120             Float_t qcl;
121             Float_t dx;
122             Float_t dz;
123           } ntuple1_st;
124
125
126           struct {
127             //      Int_t lay;
128             Int_t nx;
129             Int_t nz;
130           } ntuple2_st;
131
132
133           ntuple = new TTree("ntuple","Demo ntuple");
134           ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
135           ntuple->Branch("nx",&ntuple_st.nx,"nx/I");
136           ntuple->Branch("nz",&ntuple_st.nz,"nz/I");
137           ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
138           ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
139           ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
140           ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
141           ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
142
143           ntuple1 = new TTree("ntuple1","Demo ntuple1");
144           ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
145           ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
146           ntuple1->Branch("det",&ntuple1_st.det,"det/I");
147           ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I");
148           ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I");
149           ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F");
150           ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
151           ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
152           ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
153           ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
154
155
156           ntuple2 = new TTree("ntuple2","Demo ntuple2");
157           //      ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
158           ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I");
159           ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I");
160
161 // ------------------------------------------------------------------------
162
163           // Module loop
164
165           for (Int_t mod=0; mod<nent; mod++) {
166               AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
167
168               Int_t nhits = itsModule->GetNhits();
169               if(nhits) printf("module nhits %d %d\n",mod,nhits);
170               if(!nhits) continue;
171      
172               ITS->ResetClusters();
173               TC->GetEvent(mod);
174               Int_t nclust = ITSclusters->GetEntries();
175               if (nclust) printf("Found %d clust for module %d in det type %d \n",nclust,mod,idettype);
176               if (!nclust) continue;
177
178               // cluster/hit loops
179
180         for (Int_t clu=0;clu<nclust;clu++) {
181                 itsclu   = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu);
182                 printf("%d %d %f %f %f\n",itsclu->NclZ(),itsclu->NclX(),itsclu->Q(),itsclu->X(),itsclu->Z());
183
184                 Int_t noverlaps = 0;
185                 Int_t noverprim = 0;
186
187                 Int_t clustersizex = itsclu->NclX();
188                 Int_t clustersizez = itsclu->NclZ();
189                 //       Int_t xstart = itsclu->XStart();
190                 //       Int_t xstop = itsclu->XStop();
191                 Int_t xstart = itsclu->XStartf();
192                 Int_t xstop = itsclu->XStopf();
193                 Float_t fxstart = xstart*50;
194                 Float_t fxstop = (xstop+1)*50;
195                 Float_t zstart = itsclu->ZStart();
196                 Float_t zstop = itsclu->ZStop();
197                 Int_t zend = itsclu->Zend();
198                 Float_t clusterx = itsclu->X();
199                 Float_t clusterz = itsclu->Z();
200                 Float_t clusterQ = itsclu->Q();
201
202                 ntuple2_st.nx = clustersizex;
203                 ntuple2_st.nz = clustersizez;
204
205                 ntuple2->Fill();
206
207                 Int_t icl = 0;
208                 Float_t dxprimlast = 10.e+6;
209                 Float_t dzprimlast = 10.e+6;
210
211
212                 //        if(module > 217 && module <  226) {
213                 cout<<"mod,nclust,clu,Nxpix,Nzpix ="<<mod<<","<<nclust<<","<<clu<<","<<clustersizex<<","<<clustersizez<<endl;
214                 cout<<"clusx,clusz ="<<clusterx<<","<<clusterz<<endl;
215                 cout<<"XStartf,XStopf,ZStart,ZStop ="<<fxstart<<","<<fxstop<<","<<zstart<<","<<zstop<<endl;
216                 //         }
217                 
218
219                 Float_t SPDlength = 83600;      
220                 Float_t SPDwidth = 12800;       
221                 Float_t xhit0 = 1e+5;
222                 Float_t zhit0 = 1e+5;
223
224                 
225        for (Int_t hit=0;hit<nhits;hit++) {
226
227                   // Find coordinate differences between the hit and cluster positions
228                   // for the resolution determination.
229
230                   itsHit   = (AliITShit*)itsModule->GetHit(hit);
231
232                   Int_t hitlayer = itsHit->GetLayer();
233                   Int_t hitladder= itsHit->GetLadder();
234                   Int_t hitdet= itsHit->GetDetector();
235
236                   Int_t clusterlayer = hitlayer;
237                   Int_t clusterladder= hitladder;
238                   Int_t clusterdetector = hitdet;
239
240                   Int_t track = itsHit->fTrack;
241                   Int_t dray = 0;
242                   Int_t hitstat = itsHit->GetTrackStatus();
243
244
245                   Float_t zhit = 10000*itsHit->GetZL();
246                   Float_t xhit = 10000*itsHit->GetXL();
247
248                   if(abs(zhit) > SPDlength/2) {
249                     if(hitstat == 66) zhit0 = 1e+5;
250                     continue;
251                   }
252
253                   if(abs(xhit) > SPDwidth/2) {
254                     if(hitstat == 66) xhit0 = 1e+5;
255                     continue;
256                   }
257
258                   zhit += SPDlength/2;
259                   xhit += SPDwidth/2;
260                   Float_t yhit = 10000*itsHit->GetYL();
261
262                   if(hitlayer == 1 && hitstat == 66 && yhit > 71) {
263                     xhit0 = xhit;
264                     zhit0 = zhit;
265                   }
266                   if(hitlayer == 2 && hitstat == 66 && yhit < -71) {
267                     xhit0 = xhit;
268                     zhit0 = zhit;
269                   }
270
271                   if(hitstat != 68) continue; // Take only the hit if the last
272                   // track point went out from
273                   // the detector.
274
275                   if(xhit0 > 9e+4 || zhit0 > 9e+4) continue;
276
277
278                   Float_t xmed = (xhit + xhit0)/2;
279                   Float_t zmed = (zhit + zhit0)/2;
280
281                   Float_t xdif = xmed - clusterx;
282                   Float_t zdif = zmed - clusterz;
283
284       cout<<"clu,hit,xmed,fxstart,fxstop,zmed,zstart,zstop ="<<clu<<","<<hit<<","<<xmed<<","<<fxstart<<","<<fxstop<<","<<zmed<<","<<zstart<<","<<zstop<<endl;
285
286         // Consider the hits inside of cluster region only
287
288   if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
289
290         icl = 1;
291
292                 //        part = (TParticle *)particles.UncheckedAt(track);
293                 //        Int_t partcode = part->GetPdgCode();
294                 //              Int_t primery = gAlice->GetPrimary(track);
295
296         Int_t parent = itsHit->GetParticle()->GetFirstMother();
297         Int_t partcode = itsHit->GetParticle()->GetPdgCode();
298
299 //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
300 //                      310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
301
302
303         Float_t px = itsHit->GetPXL();
304         Float_t py = itsHit->GetPYL();
305         Float_t pz = itsHit->GetPZL();
306         Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz);
307
308               cout<<"track,partcode,pmod,parent ="<<track<<","<<partcode<<","<<pmod<<","<<parent<<endl;
309
310         Int_t hitprim = 0;
311
312         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
313                                                  // at p < 6 MeV/c
314
315         if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but
316                                                  // not for delta ray which
317                                                  // also went out from the
318                                                  // detector and returned
319                                                  // again
320
321         if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
322                                               // particles
323
324         if(hitprim > 0) noverprim = noverprim + 1;
325
326         if(hitprim > 0) {
327          dxprimlast = xdif;
328          dzprimlast = zdif;
329         }
330
331         // fill ntuple
332
333          ntuple_st.lay = hitlayer;
334          ntuple_st.nx = clustersizex;
335          ntuple_st.nz = clustersizez;
336          ntuple_st.hitprim = hitprim;
337          ntuple_st.partcode = partcode;
338          ntuple_st.dx = xdif;
339          ntuple_st.dz = zdif;
340          ntuple_st.pmod = pmod;
341
342          ntuple->Fill();
343
344       
345
346       if(hitprim > 0) {   // for primary particles
347         if(hitlayer == 1) {
348      cout<<"!!!!!! lay,hitprim,xdif,zdif ="<<hitlayer<<","<<hitprim<<","<<xdif<<","<<zdif<<endl;
349            Xres1->Fill(xdif);
350            Zres1->Fill(zdif);
351         }
352         if(hitlayer == 2) {
353      cout<<"!!!!!! lay,hitprim,xdif,zdif ="<<hitlayer<<","<<hitprim<<","<<xdif<<","<<zdif<<endl;
354            Xres2->Fill(xdif);
355            Zres2->Fill(zdif);
356         }
357       } // primery particles
358
359      } // end of cluster region
360    } // end of hit loop
361
362                 
363       if(icl == 1) {
364
365         // fill ntuple1
366
367         //      ntuple1->Fill(clusterlayer,clustersizex,clustersizez,noverlaps,\
368 noverprim,dx,dz);
369
370         if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
371                                           // delta rays only
372
373          ntuple1_st.lay = clusterlayer;
374          ntuple1_st.lad = clusterladder;
375          ntuple1_st.det = clusterdetector;
376          ntuple1_st.nx = clustersizex;
377          ntuple1_st.nz = clustersizez;
378          ntuple1_st.qcl = clusterQ;
379          ntuple1_st.noverlaps = noverlaps;
380          ntuple1_st.noverprim = noverprim;
381          ntuple1_st.dx = dxprimlast;
382          ntuple1_st.dz = dzprimlast;
383
384          ntuple1->Fill();
385
386      } // icl = 1
387     } // cluster loop
388    } // module loop       
389   } // idettype loop
390  } // end if ITS
391 } // event loop 
392
393
394    //  Write and Draw Histogramms and ntuples
395
396
397
398    TFile fhistos("SPD_his.root","RECREATE");
399
400    ntuple->Write();
401    ntuple1->Write();
402    ntuple2->Write();
403
404    Nxpix1->Write();
405    Nzpix1->Write();
406    Nxpix2->Write();
407    Nzpix2->Write();
408
409    Xpix1->Write();
410    Zpix1->Write();
411    Xpix2->Write();
412    Zpix2->Write();
413
414    Xres1->Write();
415    Zres1->Write();
416    Xres2->Write();
417    Zres2->Write();
418
419    fhistos.Close();
420    cout<<"!!! Histogramms and ntuples were written"<<endl;
421
422    TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
423    c1->Divide(2,2);
424    c1->cd(1);
425    gPad->SetFillColor(33);
426          Xres1->SetFillColor(42);
427          Xres1->Draw();
428    c1->cd(2);
429    gPad->SetFillColor(33);
430          Zres1->SetFillColor(46);
431          Zres1->Draw();
432    c1->cd(3);
433    gPad->SetFillColor(33);
434          Xres2->SetFillColor(42);
435          Xres2->Draw();
436    c1->cd(4);
437    gPad->SetFillColor(33);
438          Zres2->SetFillColor(46);
439          Zres2->Draw();
440
441 //     cout<<"END  test for clusters and hits "<<endl;
442
443 //     file->Close();   
444 }
445
446
447