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