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