]>
Commit | Line | Data |
---|---|---|
d1f98e87 | 1 | void SPDclusterTest(Int_t evNumber1=0,Int_t evNumber2=0){ |
2 | // | |
3 | // macro to monitor the SPD digitization and clusterization done with | |
4 | // the Bari/Salerno model | |
5 | // | |
6 | // R. Caliandro 15/05/2001 | |
7 | // | |
8 | // | |
9 | //--plots displayed: | |
10 | // | |
11 | //--pag1: number of hits per SPD detector (1-->250) | |
12 | // number of hits per SPD detector (1-->250) | |
13 | // number of clusters per SPD detector (1-->250) | |
14 | // | |
15 | //--pag2: r-phi cluster length layer 1 (red) | |
16 | // z cluster length layer 1 (red) | |
17 | // r-phi cluster length layer 2 (blue) | |
18 | // z cluster length layer 2 (blue) | |
19 | // | |
20 | //--pag3: r-phi resolution layer 1 (red) | |
21 | // z resolution layer 1 (red) | |
22 | // r-phi resolution layer 2 (blue) | |
23 | // z resolution layer 2 (blue) | |
24 | // | |
25 | //--pag4: Cluster shape analysis for clusters of 1, 2 and 3 digits | |
26 | // zdim versus xdim for clusters of 4 digits | |
27 | // | |
28 | // input file name, digitized and clusterized | |
29 | char *filein="galice.root"; | |
30 | // output file name, containing histograms | |
31 | char *fileout="SPD_his.root"; | |
32 | // flag for debugging: 0=no debugging, 1=debugging | |
33 | Int_t debug=0; | |
34 | ||
35 | ||
36 | // Dynamically link some shared libs | |
37 | if (gClassTable->GetID("AliRun") < 0) { | |
38 | gROOT->LoadMacro("loadlibs.C"); | |
39 | loadlibs(); | |
40 | } else { | |
41 | delete gAlice; | |
42 | gAlice=0; | |
43 | } | |
44 | ||
45 | ||
46 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
47 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filein); | |
48 | if (!file) file = new TFile(filein); | |
49 | ||
50 | // Get AliRun object from file or create it if not on file | |
51 | if (!gAlice) { | |
52 | gAlice = (AliRun*)file->Get("gAlice"); | |
53 | if (gAlice) printf("AliRun object found on file\n"); | |
54 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
55 | } | |
56 | ||
d1f98e87 | 57 | |
58 | //======================================================= | |
59 | //--booking of ntuples | |
60 | ||
61 | //--ntuple for each detector | |
62 | TNtuple *ntuple2 = new TNtuple("ntuple2","","ndet:lay:lad:det:nhits:ndig:nclus"); | |
63 | //--ntuple for each cluster | |
64 | TNtuple *ntuple = new TNtuple("ntuple","","ndet:iclus:ndigclus:xdim:zdim:xdiff:zdiff:anglex:anglez:pmom:errx:errz"); | |
65 | ||
66 | //--booking of histograms | |
67 | //layer 1 | |
68 | TH1F *hist1n1 = new TH1F("hist1n1","xdim",15,0.5,15.5); | |
69 | TH1F *hist2n1 = new TH1F("hist2n1","zdim",10,0.5,10.5); | |
70 | TH1F *hist3n1 = new TH1F("hist3n1","dig/clus",20,0.5,20.5); | |
71 | TH1F *hist4n1 = new TH1F("hist4n1","errx",100,0,0.01); | |
72 | TH1F *hist5n1 = new TH1F("hist5n1","errz",500,0,0.05); | |
73 | TH2F *hist7n1 = new TH2F("hist7n1","xdim:delx",80,0,800.,15,0.5,15.5); | |
74 | TH2F *hist8n1 = new TH2F("hist8n1","zdim:delz",180,0,1800.,10,0.5,10.5); | |
75 | //layer 2 | |
76 | TH1F *hist1n2 = new TH1F("hist1n2","xdim",15,0.5,15.5); | |
77 | TH1F *hist2n2 = new TH1F("hist2n2","zdim",10,0.5,10.5); | |
78 | TH1F *hist3n2 = new TH1F("hist3n2","dig/clus",20,0.5,20.5); | |
79 | TH1F *hist4n2 = new TH1F("hist4n2","errx",100,0,0.01); | |
80 | TH1F *hist5n2 = new TH1F("hist5n2","errz",500,0,0.05); | |
81 | TH2F *hist7n2 = new TH2F("hist7n2","xdim:delx",80,0,800.,15,0.5,15.5); | |
82 | TH2F *hist8n2 = new TH2F("hist8n2","zdim:delz",180,0,1800.,10,0.5,10.5); | |
83 | //--resolution | |
84 | TH1F *hist1 = new TH1F("hist1","xdiff",200,-100,100); | |
85 | TH1F *hist3 = new TH1F("hist3","xdiff",200,-100,100); | |
86 | TH1F *hist2 = new TH1F("hist2","zdiff",170,-850,850); | |
87 | TH1F *hist4 = new TH1F("hist4","zdiff",170,-850,850); | |
88 | //--momentum | |
89 | TH1F *hist5 = new TH1F("hist5","pmom",200,0,2000); | |
90 | //--rapidity | |
91 | TH1F *hist6 = new TH1F("hist6","rapidity",60,-3,3); | |
92 | TH1F *hist6b= new TH1F("hist6b","rapidity - charged tracks",60,-3,3); | |
93 | TH1F *hist6b1= new TH1F("hist6b1","rapidity - charged tracks SPD",60,-3,3); | |
94 | //--pseudo-rapidity | |
95 | TH1F *hist6p = new TH1F("hist6p","eta - charged tracks ",60,-3,3); | |
96 | TH1F *hist6p1 = new TH1F("hist6p1","eta - charged tracks SPD ",60,-3,3); | |
97 | TH1F *hist6p2 = new TH1F("hist6p2","eta - charged tracks SPD 2 ",60,-3,3); | |
98 | //--resolution vs angle | |
99 | TH1F *hist11n1=new TH1F("hist11n1","anglex - layer 1",180,-90,90); | |
100 | TH1F *hist11n2=new TH1F("hist11n2","anglex - layer 2",180,-90,90); | |
101 | TH1F *hist12n1=new TH1F("hist12n1","anglez - layer 1",360,-180,180); | |
102 | TH1F *hist12n2=new TH1F("hist12n2","anglez - layer 2",360,-180,180); | |
103 | TH2F *hist13n1=new TH2F("hist13n1","xidff:anglex",20,-15,15,200,-100,100); | |
104 | TH2F *hist13n2=new TH2F("hist13n2","xidff:anglex",20,-30,-15,200,-100,100); | |
105 | TH2F *hist14n1=new TH2F("hist14n1","zidff:anglez",18,-90,90,170,-850,850); | |
106 | TH2F *hist14n2=new TH2F("hist14n2","zidff:anglez",18,-90,90,170,-850,850); | |
107 | //--histograms for cluster shape analysis | |
108 | TH1F *histsp1=new TH1F("histsp1","Cluster shape (1)",10,0.5,10.5); | |
109 | TH2F *histsp2=new TH2F("histsp2","Cluster shape (2)",5,0.5,5.5,5,0.5,5.5); | |
110 | //======================================================= | |
111 | ||
112 | //loop over events | |
113 | for (int nev=0; nev<= evNumber2; nev++) { | |
114 | Int_t nparticles = gAlice->GetEvent(nev); | |
115 | cout << "nev " <<nev<<endl; | |
116 | cout << "nparticles " <<nparticles<<endl; | |
117 | if (nev < evNumber1) continue; | |
118 | if (nparticles <= 0) return; | |
119 | ||
120 | TTree *TH = gAlice->TreeH(); | |
121 | Int_t ntracks = TH->GetEntries(); | |
122 | cout << "ntracks " <<ntracks<<endl; | |
123 | ||
124 | // Get pointers to Alice detectors and Digit containers | |
125 | AliITS *ITS = (AliITS *)gAlice->GetModule("ITS"); | |
4d655c07 | 126 | //to get the segmentation pointer |
127 | AliITSDetType *iDetType=ITS->DetType(0); | |
128 | AliITSsegmentationSPD *seg=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); | |
129 | //TClonesArray *Particles = gAlice->Particles(); | |
d1f98e87 | 130 | if(!ITS) return; |
131 | ||
132 | // fill modules with sorted by module hits | |
133 | Int_t nmodules; | |
134 | ITS->InitModules(-1,nmodules); | |
135 | ITS->FillModules(nev,evNumber2,nmodules," "," "); | |
136 | ||
137 | // get pointer to modules array | |
138 | TObjArray *mods = ITS->GetModules(); | |
139 | AliITShit *itsHit; | |
140 | ||
141 | ||
142 | //get the Tree for clusters | |
143 | ITS->GetTreeC(nev); | |
144 | TTree *TC=ITS->TreeC(); | |
145 | //TC->Print(); | |
146 | Int_t nent=TC->GetEntries(); | |
147 | printf("Found %d entries in the tree of clusters)\n",nent); | |
148 | TClonesArray *ITSclusters = ITS->ClustersAddress(0); | |
149 | printf("ITSclusters %p\n",ITSclusters); | |
150 | ||
151 | //get the Tree for digits | |
152 | TTree *TD = gAlice->TreeD(); | |
153 | //TD->Print(); | |
154 | Int_t nentd=TD->GetEntries(); | |
155 | printf("Found %d entries in the tree of digits)\n",nentd); | |
156 | TObjArray *fBranches=TD->GetListOfBranches(); | |
157 | TBranch *branch = (TBranch*)fBranches->UncheckedAt(0); | |
158 | printf ("branch %p entries %d \n",branch,branch->GetEntries()); | |
159 | TClonesArray *ITSdigits = ITS->DigitsAddress(0); | |
160 | printf ("ITSdigits %p \n",ITSdigits); | |
161 | ||
162 | //get the Tree for rec points | |
163 | TTree *TR = gAlice->TreeR(); | |
164 | //TR->Print(); | |
165 | Int_t nentr=TR->GetEntries(); | |
166 | printf("Found %d entries in the tree of rec points)\n",nentr); | |
167 | TClonesArray *ITSrec = ITS->RecPoints(); | |
168 | printf ("ITSrec %p \n",ITSrec); | |
169 | AliITSRecPoint *recp; | |
170 | ||
171 | // calculus of rapidity distribution for the generated tracks | |
172 | gAlice-> ResetHits(); | |
173 | TParticle *particle; | |
174 | for (Int_t track=0; track<ntracks; track++) | |
175 | { | |
176 | particle = (TParticle*)gAlice->Particle(track); | |
177 | Int_t ikparen = particle -> GetFirstMother(); | |
178 | Double_t charge = particle -> GetPDG() ->Charge(); | |
179 | charge = charge/3.; //charge is multiplied by 3 in PDG | |
180 | Double_t mass = particle -> GetPDG() -> Mass(); | |
181 | Double_t eta = particle -> Eta(); | |
182 | Int_t pdgcode = particle -> GetPdgCode(); | |
183 | char* title = particle -> GetTitle(); | |
184 | if (ikparen<0) | |
185 | { | |
186 | Double_t part_ene = particle->Energy(); | |
187 | Double_t part_pz = particle->Pz(); | |
188 | Double_t rapid; | |
189 | if (part_ene != part_pz) | |
190 | { | |
191 | rapid=0.5*TMath::Log((part_ene+part_pz)/(part_ene-part_pz)); | |
192 | } | |
193 | else { | |
194 | rapid = 1.e30; | |
195 | } | |
196 | // filling of the rapidity histogram | |
197 | hist6->Fill( (Float_t) rapid); | |
198 | if( charge != 0 ) { | |
199 | hist6b->Fill( (Float_t) rapid); | |
200 | hist6p->Fill( (Float_t) eta); | |
201 | } | |
202 | // printf("charge= %f, mass = %f , pdg= %d, title = %s\n", | |
203 | // charge,mass,pdgcode,title); | |
204 | } | |
205 | } | |
206 | ||
207 | AliITSgeom *g = ((AliITS *)ITS)->GetITSgeom(); | |
208 | Int_t lay, lad, det; | |
209 | //printf("Starts loop on SPD detectors\n"); | |
210 | ||
211 | ||
212 | //loop over the pixel detectors index=0-79 (1-20)*4 layer 1 | |
213 | // index=80-239 (1-40)*4 layer 2 | |
214 | for (Int_t index=g->GetStartSPD();index<=g->GetLastSPD();index++) | |
215 | // for (Int_t index=g->GetStartSPD();index<1;index++) //debug | |
216 | { | |
217 | ||
218 | g->GetModuleId(index,lay,lad,det); | |
219 | //printf("detector %d (lay=%d lad=%d det=%d)\n",index+1,lay,lad,det); | |
220 | ||
221 | AliITSmodule *itsModule = (AliITSmodule*) mods->At(index); | |
222 | Int_t numofhits = itsModule->GetNhits(); | |
223 | //printf("number of hits %d\n",numofhits); | |
224 | if(!numofhits) continue; | |
225 | ||
226 | //---------- starts test on digits | |
227 | ITS->ResetDigits(); | |
228 | TD->GetEvent(index); | |
229 | Int_t ndigits = ITSdigits->GetEntriesFast(); | |
230 | //if (ndigits) printf("Found %d digits for module %d \n",ndigits,index+1); | |
231 | if (!ndigits) printf("no digits found \n"); | |
232 | ||
233 | ||
234 | ||
235 | if(debug==1) { | |
236 | //loop on digits | |
237 | for (Int_t digit=0;digit<ndigits;digit++) { | |
238 | ITSdigit = (AliITSdigitSPD*)ITSdigits->UncheckedAt(digit); | |
239 | printf("digit=%d fCoord1=%d FCoord2=%d fSignal=%d fTracks=%d fHits=%d \n",digit,ITSdigit->fCoord1,ITSdigit->fCoord2,ITSdigit->fSignal,ITSdigit->fTracks[0],ITSdigit->fHits[0]); | |
240 | } | |
241 | cout<<"END test for digits "<<endl; | |
242 | } | |
243 | ||
244 | ||
245 | //---------- starts test on clusters | |
246 | ITS->ResetClusters(); | |
247 | TC->GetEvent(index); | |
248 | Int_t nclust = ITSclusters->GetEntries(); | |
249 | //printf("Found %d clusters \n",nclust); | |
250 | if (!nclust) printf("no clusters found \n"); | |
251 | ||
252 | ||
253 | if(debug==1) { | |
254 | //loop on clusters | |
255 | for (Int_t clu=0;clu<nclust;clu++) | |
256 | { | |
257 | //itsclu = (AliITSRawClusterSPD*) ITSclusters->UncheckedAt(clu); | |
258 | itsclu = (AliITSRawClusterSPD*) ITSclusters->At(clu); | |
259 | printf("cluster %d nZ=%f nX=%f Q=%f Z=%f X=%f\n",clu+1,itsclu->NclZ(), | |
260 | itsclu->NclX(),itsclu->Q(),itsclu->Z(),itsclu->X()); | |
261 | } | |
262 | cout<<"END test for clusters "<<endl; | |
263 | } | |
264 | ||
265 | ||
266 | ||
267 | //---------- starts test on rec points | |
268 | ITS->ResetRecPoints(); | |
269 | TR->GetEvent(index); | |
270 | Int_t nrecpoints = ITSrec->GetEntries(); | |
271 | //printf("Found %d recpoints for module %d \n",nrecpoints,index+1); | |
272 | if (!nrecpoints) printf("no recpoints found \n"); | |
273 | ||
274 | if(debug==1) { | |
275 | //loop on rec points | |
276 | for (Int_t irec=0;irec<nrecpoints;irec++) { | |
277 | recp = (AliITSRecPoint*)ITSrec->UncheckedAt(irec); | |
278 | printf("%d %f %f %f %f %d %d %d\n",irec+1,recp->GetX(),recp->GetZ(), | |
279 | recp->fSigmaX2,recp->fSigmaZ2, | |
280 | recp->fTracks[0],recp->fTracks[1],recp->fTracks[2]); | |
281 | } | |
282 | } | |
283 | ||
284 | printf("Detector n.%d (%d hits) (%d digits) (%d clusters)\n", | |
285 | index+1,numofhits,ndigits,nclust); | |
286 | ||
287 | // fill ntuple2 | |
288 | ntuple2->Fill ( (Float_t) index+1, | |
289 | (Float_t) lay, | |
290 | (Float_t) lad, | |
291 | (Float_t) det, | |
292 | (Float_t) numofhits, | |
293 | (Float_t) ndigits, | |
294 | (Float_t) nclust); | |
295 | ||
296 | Int_t xlow; | |
297 | Int_t zlow; | |
298 | Int_t xhigh; | |
299 | Int_t zhigh; | |
300 | Int_t colcenter; | |
301 | Int_t rowcenter; | |
302 | ||
303 | // loop on clusters in each detector | |
304 | for (Int_t i=0; i<nclust; i++) | |
305 | { | |
306 | ||
307 | irawclu = (AliITSRawClusterSPD*) ITSclusters->UncheckedAt(i); | |
308 | irecp = (AliITSRecPoint*)ITSrec->UncheckedAt(i); | |
309 | ||
310 | Int_t xdim = irawclu->NclX(); | |
311 | Int_t zdim = irawclu->NclZ(); | |
312 | Float_t errx = TMath::Sqrt(irecp->fSigmaX2); | |
313 | Float_t errz = TMath::Sqrt(irecp->fSigmaZ2); | |
314 | Float_t xcenter = irawclu->X(); | |
315 | Float_t zcenter = irawclu->Z(); | |
316 | Float_t ndigclus = irawclu->Q(); | |
317 | Int_t itrackclus = irecp->fTracks[0]; | |
318 | ||
319 | //Find the hits associated to the main track of the cluster | |
320 | // loop on hits in the detector | |
321 | for (Int_t hit=0; hit<numofhits; hit++) | |
322 | { | |
323 | AliITShit *itsHit = (AliITShit*)itsModule->GetHit(hit); | |
324 | Int_t itrackhit = itsHit->GetTrack(); | |
325 | //Take the same track index of the main track of the cluster | |
326 | if (itrackhit == itrackclus) { | |
327 | if (itsHit->GetTrackStatus()==66) { | |
328 | Float_t x1l = 10000*itsHit->GetXL(); //in microns | |
329 | Float_t y1l = 10000*itsHit->GetYL(); | |
330 | Float_t z1l = 10000*itsHit->GetZL(); | |
331 | Float_t p1x = 1000*itsHit->GetPXL(); //in MeV/c | |
332 | Float_t p1y = 1000*itsHit->GetPYL(); | |
333 | Float_t p1z = 1000*itsHit->GetPZL(); | |
334 | } | |
335 | else { | |
336 | Float_t x2l = 10000*itsHit->GetXL(); | |
337 | Float_t y2l = 10000*itsHit->GetYL(); | |
338 | Float_t z2l = 10000*itsHit->GetZL(); | |
339 | Float_t p2x = 1000*itsHit->GetPXL(); | |
340 | Float_t p2y = 1000*itsHit->GetPYL(); | |
341 | Float_t p2z = 1000*itsHit->GetPZL(); | |
342 | ||
343 | ||
344 | } | |
345 | } | |
346 | }// end loop on hits on detector | |
347 | ||
348 | ||
349 | Float_t pmom=TMath::Sqrt(p1x*p1x+p1y*p1y+p1z*p1z); | |
350 | hist5->Fill(pmom); | |
351 | ||
352 | Float_t dxhit = TMath::Abs(x2l-x1l); | |
353 | Float_t dzhit = TMath::Abs(z2l-z1l); | |
354 | ||
355 | Float_t xmidhit = (x1l + x2l)/2; | |
356 | Float_t zmidhit = (z1l + z2l)/2; | |
357 | ||
358 | // printf("cluster n.%d: x=%f z=%f\n",i,xcenter,zcenter); | |
359 | // printf("track n.%d: x1=%f x2=%f z1=%f z2=%f\n",itrackclus, | |
360 | // x1l, x2l, z1l, z2l); | |
361 | ||
362 | // analysis of resolution vs angle | |
363 | if(index<80) | |
364 | { | |
365 | Float_t px = -p1x; | |
366 | Float_t py = -p1y; | |
367 | } | |
368 | else{ | |
369 | Float_t px = p1x; | |
370 | Float_t py = p1y; | |
371 | } | |
372 | Float_t pz = p1z; | |
373 | // anglex is the angle in xy plane (local frame) | |
374 | Float_t anglex = atan2(px,py); | |
375 | // anglez is the angle in zy plane (local frame) | |
376 | Float_t anglez = atan2(pz,py); | |
377 | anglex *= 180.0/TMath::Pi(); // degrees | |
378 | anglez *= 180.0/TMath::Pi(); // degrees | |
379 | ||
380 | if(xmidhit != 0 || zmidhit != 0) | |
381 | { | |
382 | Float_t xdiff = (xcenter - xmidhit); | |
383 | Float_t zdiff = (zcenter - zmidhit); | |
384 | ||
385 | if(index<80) | |
386 | { | |
387 | // resolution plots | |
388 | hist1->Fill(xdiff); | |
389 | hist2->Fill(zdiff); | |
390 | ||
391 | // plots of resolution vs angle | |
392 | hist11n1->Fill(anglex); | |
393 | hist12n1->Fill(anglez); | |
394 | hist13n1->Fill(anglex,xdiff); | |
395 | hist14n1->Fill(anglez,zdiff); | |
396 | ||
397 | } else { | |
398 | ||
399 | // resolution plots | |
400 | hist3->Fill(xdiff); | |
401 | hist4->Fill(zdiff); | |
402 | ||
403 | // plots of resolution vs angle | |
404 | hist11n2->Fill(anglex); | |
405 | hist12n2->Fill(anglez); | |
406 | hist13n2->Fill(anglex,xdiff); | |
407 | hist14n2->Fill(anglez,zdiff); | |
408 | ||
409 | } | |
410 | } | |
411 | ||
412 | // fill the ntuple | |
413 | ntuple->Fill ( (Float_t) index, | |
414 | (Float_t) i, | |
415 | (Float_t) ndigclus, | |
416 | (Float_t) xdim, | |
417 | (Float_t) zdim, | |
418 | (Float_t) xdiff, | |
419 | (Float_t) zdiff, | |
420 | (Float_t) anglex, | |
421 | (Float_t) anglez, | |
422 | (Float_t) pmom, | |
423 | errx, | |
424 | errz); | |
425 | ||
426 | // other histograms | |
427 | if(index<80) | |
428 | { | |
429 | hist1n1->Fill((Float_t) xdim); | |
430 | hist2n1->Fill((Float_t) zdim); | |
431 | hist3n1->Fill(ndigclus); | |
432 | hist4n1->Fill(errx); | |
433 | hist5n1->Fill(errz); | |
434 | hist7n1->Fill(dxhit,(Float_t) xdim); | |
435 | hist8n1->Fill(dzhit,(Float_t) zdim); | |
436 | ||
437 | } else { | |
438 | hist1n2->Fill((Float_t) xdim); | |
439 | hist2n2->Fill((Float_t) zdim); | |
440 | hist3n2->Fill(ndigclus); | |
441 | hist4n2->Fill(errx); | |
442 | hist5n2->Fill(errz); | |
443 | hist7n2->Fill(dxhit,(Float_t) xdim); | |
444 | hist8n2->Fill(dzhit,(Float_t) zdim); | |
445 | } | |
446 | ||
447 | //histograms for cluster shape analysis | |
448 | Int_t xx; | |
449 | if(ndigclus<=3) { | |
450 | if(ndigclus==1) { | |
451 | xx=1; | |
452 | } elseif(ndigclus==2){ | |
453 | if(zdim==2 && xdim==1) xx=2; | |
454 | if(zdim==1 && xdim==2) xx=3; | |
455 | if(zdim==2 && xdim==2) xx=4; | |
456 | } elseif(ndigclus==3){ | |
457 | if(zdim==2 && xdim==2) xx=5; | |
458 | if(zdim==3 && xdim==1) xx=6; | |
459 | if(zdim==1 && xdim==3) xx=7; | |
460 | if(zdim==3 && xdim==3) xx=8; | |
461 | if(zdim==3 && xdim==2) xx=9; | |
462 | if(zdim==2 && xdim==3) xx=10; | |
463 | } | |
464 | histsp1->Fill((Float_t) xx); | |
465 | } elseif(ndigclus==4){ | |
466 | histsp2->Fill((Float_t) xdim,(Float_t) zdim); | |
467 | } | |
468 | ||
469 | ||
470 | }//end loop on clusters | |
471 | ||
472 | } //end loop on the SPD detectors | |
473 | ||
474 | } //end loop on events | |
475 | ||
476 | ||
477 | //================== Plot the results =================================== | |
478 | ||
479 | gROOT->Reset(); | |
480 | gStyle->SetFillColor(0); | |
481 | gStyle->SetStatW(0.37); | |
482 | gStyle->SetStatH(0.22); | |
483 | ||
484 | //----------------------------------------------------- page 1 | |
485 | gStyle->SetOptLogy(0); | |
486 | gStyle->SetOptStat(1100); | |
487 | ||
488 | TCanvas *c1 = new TCanvas("c1","hits, digits, clusters",200,10,700,780); | |
489 | c1->SetFillColor(0); | |
490 | c1->Divide(1,3); | |
491 | c1->cd(1); | |
492 | ntuple2->SetMarkerColor(kRed); | |
493 | ntuple2->SetMarkerStyle(20); | |
494 | ntuple2->SetMarkerSize(0.6); | |
495 | ntuple2->Draw("nhits:ndet",""); | |
496 | c1->cd(2); | |
497 | ntuple2->SetMarkerColor(kBlue); | |
498 | ntuple2->Draw("ndig:ndet",""); | |
499 | c1->cd(3); | |
500 | ntuple2->SetMarkerColor(6); | |
501 | ntuple2->Draw("nclus:ndet",""); | |
502 | ||
503 | //----------------------------------------------------- page 2 | |
504 | ||
505 | TCanvas *c2 = new TCanvas("c2","Cluster Lengths",200,10,700,780); | |
506 | // | |
507 | // Inside this canvas, we create 4 pads | |
508 | // | |
509 | pad1 = new TPad("pad1","xdim layer1" ,0.01,0.51,0.49,0.99,10); | |
510 | pad2 = new TPad("pad2","zdim layer1" ,0.51,0.51,0.99,0.99,10); | |
511 | pad3 = new TPad("pad3","xdim layer2" ,0.01,0.01,0.49,0.49,10); | |
512 | pad4 = new TPad("pad4","zdim layer2" ,0.51,0.01,0.99,0.49,10); | |
513 | pad1->Draw(); | |
514 | pad2->Draw(); | |
515 | pad3->Draw(); | |
516 | pad4->Draw(); | |
517 | // | |
518 | gStyle->SetStatW(0.40); | |
519 | gStyle->SetStatH(0.20); | |
520 | gStyle->SetStatColor(42); | |
521 | gStyle->SetOptStat(111110); | |
522 | // gStyle->SetOptFit(1); | |
523 | ||
524 | pad1->cd(); | |
525 | pad1->SetLogy(); | |
526 | hist1n1->Draw(); | |
527 | hist1n1->SetLineWidth(2); | |
528 | hist1n1->SetLineColor(kRed); | |
529 | hist1n1->GetXaxis()->SetNdivisions(110); | |
530 | hist1n1->GetYaxis()->SetLabelSize(0.06); | |
531 | c2->Update(); | |
532 | // | |
533 | pad2->cd(); | |
534 | pad2->SetLogy(); | |
535 | hist2n1->Draw(); | |
536 | hist2n1->SetLineWidth(2); | |
537 | hist2n1->SetLineColor(kRed); | |
538 | hist2n1->GetXaxis()->SetNdivisions(110); | |
539 | hist2n1->GetYaxis()->SetLabelSize(0.06); | |
540 | c2->Update(); | |
541 | // | |
542 | pad3->cd(); | |
543 | pad3->SetLogy(); | |
544 | hist1n2->Draw(); | |
545 | hist1n2->SetLineWidth(2); | |
546 | hist1n2->SetLineColor(kBlue); | |
547 | hist1n2->GetXaxis()->SetNdivisions(110); | |
548 | hist1n2->GetYaxis()->SetLabelSize(0.06); | |
549 | c2->Update(); | |
550 | // | |
551 | pad4->cd(); | |
552 | pad4->SetLogy(); | |
553 | hist2n2->Draw(); | |
554 | hist2n2->SetLineColor(kBlue); | |
555 | hist2n2->SetLineWidth(2); | |
556 | hist2n2->GetXaxis()->SetNdivisions(110); | |
557 | hist2n2->GetYaxis()->SetLabelSize(0.06); | |
558 | c2->Update(); | |
559 | // | |
560 | //----------------------------------------------------- page 3 | |
561 | ||
562 | TCanvas *c3 = new TCanvas("c3","Resolutions",200,10,700,780); | |
563 | // | |
564 | // Inside this canvas, we create 4 pads | |
565 | // | |
566 | pad1 = new TPad("pad1","xdiff layer1" ,0.01,0.51,0.49,0.99,10); | |
567 | pad2 = new TPad("pad2","zdiff layer1" ,0.51,0.51,0.99,0.99,10); | |
568 | pad3 = new TPad("pad3","xdiff layer2" ,0.01,0.01,0.49,0.49,10); | |
569 | pad4 = new TPad("pad4","zdiff layer2" ,0.51,0.01,0.99,0.49,10); | |
570 | pad1->Draw(); | |
571 | pad2->Draw(); | |
572 | pad3->Draw(); | |
573 | pad4->Draw(); | |
574 | // | |
575 | gStyle->SetStatW(0.20); | |
576 | gStyle->SetStatH(0.20); | |
577 | gStyle->SetStatColor(42); | |
578 | gStyle->SetOptStat(0); | |
579 | gStyle->SetOptFit(1); | |
580 | ||
581 | pad1->cd(); | |
582 | hist1->Draw(); | |
583 | hist1->SetLineColor(kRed); | |
584 | hist1->Fit("gaus"); | |
585 | c3->Update(); | |
586 | // | |
587 | pad2->cd(); | |
588 | hist2->Draw(); | |
589 | hist2->SetLineColor(kRed); | |
590 | hist2->Fit("gaus"); | |
591 | c3->Update(); | |
592 | // | |
593 | pad3->cd(); | |
594 | hist3->Draw(); | |
595 | hist3->SetLineColor(kBlue); | |
596 | hist3->Fit("gaus"); | |
597 | c3->Update(); | |
598 | // | |
599 | pad4->cd(); | |
600 | hist4->Draw(); | |
601 | hist4->SetLineColor(kBlue); | |
602 | hist4->Fit("gaus"); | |
603 | c3->Update(); | |
604 | ||
605 | //----------------------------------------------------- page 4 | |
606 | TCanvas *c4 = new TCanvas("c4","Cluster Shape Analysis",200,10,700,780); | |
607 | // | |
608 | gStyle->SetOptStat(0); | |
609 | c4->SetFillColor(0); | |
610 | c4->Divide(1,2); | |
611 | ||
612 | c4->cd(1); | |
613 | c4_1->SetLogy(); | |
614 | histsp1->Draw(); | |
615 | c4->cd(2); | |
616 | c4_2->Divide(2,1); | |
617 | c4_2->cd(1); | |
618 | histsp2->Draw("box"); | |
619 | c4_2->cd(2); | |
620 | clustershape(); | |
621 | ||
622 | ||
623 | ||
624 | ||
625 | //================== Store the histograms =================================== | |
626 | ||
627 | //to write the histograms into a .root file | |
628 | TFile outfile(fileout,"RECREATE"); | |
629 | ||
630 | ntuple->Write(); | |
631 | ntuple2->Write(); | |
632 | hist1n1->Write(); | |
633 | hist2n1->Write(); | |
634 | hist3n1->Write(); | |
635 | hist4n1->Write(); | |
636 | hist5n1->Write(); | |
637 | hist7n1->Write(); | |
638 | hist8n1->Write(); | |
639 | hist1n2->Write(); | |
640 | hist2n2->Write(); | |
641 | hist3n2->Write(); | |
642 | hist4n2->Write(); | |
643 | hist5n2->Write(); | |
644 | hist7n2->Write(); | |
645 | hist8n2->Write(); | |
646 | hist1->Write(); | |
647 | hist2->Write(); | |
648 | hist3->Write(); | |
649 | hist4->Write(); | |
650 | hist5->Write(); | |
651 | hist6->Write(); | |
652 | hist6b->Write(); | |
653 | hist6p->Write(); | |
654 | hist6b1->Write(); | |
655 | hist6p1->Write(); | |
656 | hist6p2->Write(); | |
657 | hist11n1->Write(); | |
658 | hist11n2->Write(); | |
659 | hist12n1->Write(); | |
660 | hist12n2->Write(); | |
661 | hist13n1->Write(); | |
662 | hist13n2->Write(); | |
663 | hist14n1->Write(); | |
664 | hist14n2->Write(); | |
665 | histsp1->Write(); | |
666 | histsp2->Write(); | |
667 | ||
668 | outfile->Close(); | |
669 | } | |
670 | //----------------------------------------------------------------- | |
671 | ||
672 | ||
673 | ||
674 | void clustershape(){ | |
675 | // | |
676 | //macro to display the legend of the cluster shape analysis plot | |
677 | // | |
678 | ||
679 | TPad *pad1 = new TPad("pad1", "This is pad2",0,0,1,1); | |
680 | pad1->Draw(); | |
681 | pad1->cd(); | |
682 | pad1->Range(0,0.25,1,1); | |
683 | pad1->SetFillColor(0); | |
684 | pad1->SetBorderSize(1); | |
685 | ||
686 | //------------------------------------------ | |
687 | Float_t yfirst= 0.95; | |
688 | Float_t ysh = 0.05; | |
689 | Float_t ysiz = 0.02; | |
690 | Float_t ynum = 0.005; | |
691 | //------------------------------------------ | |
692 | ||
693 | //bin 1 | |
694 | TLatex *tex = new TLatex(0.12,yfirst,"1"); | |
695 | tex->SetTextSize(0.07); | |
696 | tex->SetLineWidth(2); | |
697 | tex->Draw(); | |
698 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
699 | pave->SetFillColor(18); | |
700 | pave->SetLineWidth(2); | |
701 | pave->Draw(); | |
702 | ||
703 | //bin 2 | |
704 | yfirst=yfirst-ysh; | |
705 | TLatex *tex = new TLatex(0.12,yfirst,"2"); | |
706 | tex->SetTextSize(0.07); | |
707 | tex->SetLineWidth(2); | |
708 | tex->Draw(); | |
709 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
710 | pave->SetFillColor(18); | |
711 | pave->SetLineWidth(2); | |
712 | pave->Draw(); | |
713 | TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br"); | |
714 | pave->SetFillColor(18); | |
715 | pave->SetLineWidth(2); | |
716 | pave->Draw(); | |
717 | ||
718 | //bin 3 | |
719 | yfirst=yfirst-ysh; | |
720 | TLatex *tex = new TLatex(0.12,yfirst-ynum,"3"); | |
721 | tex->SetTextSize(0.07); | |
722 | tex->SetLineWidth(2); | |
723 | tex->Draw(); | |
724 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
725 | pave->SetFillColor(18); | |
726 | pave->SetLineWidth(2); | |
727 | pave->Draw(); | |
728 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst-ysiz,1,"br"); | |
729 | pave->SetFillColor(18); | |
730 | pave->SetLineWidth(2); | |
731 | pave->Draw(); | |
732 | ||
733 | //bin 4 | |
734 | yfirst=yfirst-1.8*ysh; | |
735 | TLatex *tex = new TLatex(0.12,yfirst+3*ynum,"4"); | |
736 | tex->SetTextSize(0.07); | |
737 | tex->SetLineWidth(2); | |
738 | tex->Draw(); | |
739 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
740 | pave->SetFillColor(18); | |
741 | pave->SetLineWidth(2); | |
742 | pave->Draw(); | |
743 | TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br"); | |
744 | pave->SetFillColor(18); | |
745 | pave->SetLineWidth(2); | |
746 | pave->Draw(); | |
747 | ||
748 | //bin 5 | |
749 | yfirst=yfirst-1.5*ysh; | |
750 | TLatex *tex = new TLatex(0.12,yfirst+3*ynum,"5"); | |
751 | tex->SetTextSize(0.07); | |
752 | tex->SetLineWidth(2); | |
753 | tex->Draw(); | |
754 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
755 | pave->SetFillColor(18); | |
756 | pave->SetLineWidth(2); | |
757 | pave->Draw(); | |
758 | TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br"); | |
759 | pave->SetFillColor(18); | |
760 | pave->SetLineWidth(2); | |
761 | pave->Draw(); | |
762 | TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br"); | |
763 | pave->SetFillColor(18); | |
764 | pave->SetLineWidth(2); | |
765 | pave->Draw(); | |
766 | ||
767 | //bin 6 | |
768 | yfirst=yfirst-1.5*ysh; | |
769 | TLatex *tex = new TLatex(0.12,yfirst+ynum,"6"); | |
770 | tex->SetTextSize(0.07); | |
771 | tex->SetLineWidth(2); | |
772 | tex->Draw(); | |
773 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
774 | pave->SetFillColor(18); | |
775 | pave->SetLineWidth(2); | |
776 | pave->Draw(); | |
777 | TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br"); | |
778 | pave->SetFillColor(18); | |
779 | pave->SetLineWidth(2); | |
780 | pave->Draw(); | |
781 | TPave *pave = new TPave(0.7,yfirst,0.9,yfirst+ysiz,1,"br"); | |
782 | pave->SetFillColor(18); | |
783 | pave->SetLineWidth(2); | |
784 | pave->Draw(); | |
785 | ||
786 | //bin 7 | |
787 | yfirst=yfirst-2*ysh; | |
788 | TLatex *tex = new TLatex(0.12,yfirst+ysiz+ynum,"7"); | |
789 | tex->SetTextSize(0.07); | |
790 | tex->SetLineWidth(2); | |
791 | tex->Draw(); | |
792 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
793 | pave->SetFillColor(18); | |
794 | pave->SetLineWidth(2); | |
795 | pave->Draw(); | |
796 | TPave *pave = new TPave(0.3,yfirst+ysiz,0.5,yfirst+2*ysiz,1,"br"); | |
797 | pave->SetFillColor(18); | |
798 | pave->SetLineWidth(2); | |
799 | pave->Draw(); | |
800 | TPave *pave = new TPave(0.3,yfirst+2*ysiz,0.5,yfirst+3*ysiz,1,"br"); | |
801 | pave->SetFillColor(18); | |
802 | pave->SetLineWidth(2); | |
803 | pave->Draw(); | |
804 | ||
805 | //bin 8 | |
806 | yfirst=yfirst-1.5*ysh; | |
807 | TLatex *tex = new TLatex(0.12,yfirst+ynum,"8"); | |
808 | tex->SetTextSize(0.07); | |
809 | tex->SetLineWidth(2); | |
810 | tex->Draw(); | |
811 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
812 | pave->SetFillColor(18); | |
813 | pave->SetLineWidth(2); | |
814 | pave->Draw(); | |
815 | TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br"); | |
816 | pave->SetFillColor(18); | |
817 | pave->SetLineWidth(2); | |
818 | pave->Draw(); | |
819 | TPave *pave = new TPave(0.7,yfirst+2*ysiz,0.9,yfirst+3*ysiz,1,"br"); | |
820 | pave->SetFillColor(18); | |
821 | pave->SetLineWidth(2); | |
822 | pave->Draw(); | |
823 | ||
824 | //bin 9 | |
825 | yfirst=yfirst-1.5*ysh; | |
826 | TLatex *tex = new TLatex(0.12,yfirst+ynum,"9"); | |
827 | tex->SetTextSize(0.07); | |
828 | tex->SetLineWidth(2); | |
829 | tex->Draw(); | |
830 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
831 | pave->SetFillColor(18); | |
832 | pave->SetLineWidth(2); | |
833 | pave->Draw(); | |
834 | TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br"); | |
835 | pave->SetFillColor(18); | |
836 | pave->SetLineWidth(2); | |
837 | pave->Draw(); | |
838 | TPave *pave = new TPave(0.7,yfirst+ysiz,0.9,yfirst+2*ysiz,1,"br"); | |
839 | pave->SetFillColor(18); | |
840 | pave->SetLineWidth(2); | |
841 | pave->Draw(); | |
842 | ||
843 | //bin 10 | |
844 | yfirst=yfirst-1.5*ysh; | |
845 | TLatex *tex = new TLatex(0.12,yfirst-ynum,"10"); | |
846 | tex->SetTextSize(0.07); | |
847 | tex->SetLineWidth(2); | |
848 | tex->Draw(); | |
849 | TPave *pave = new TPave(0.3,yfirst-ysiz,0.5,yfirst,1,"br"); | |
850 | pave->SetFillColor(18); | |
851 | pave->SetLineWidth(2); | |
852 | pave->Draw(); | |
853 | TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br"); | |
854 | pave->SetFillColor(18); | |
855 | pave->SetLineWidth(2); | |
856 | pave->Draw(); | |
857 | TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br"); | |
858 | pave->SetFillColor(18); | |
859 | pave->SetLineWidth(2); | |
860 | pave->Draw(); | |
861 | } | |
862 |