Updated version of the Bari code to work with the HEAD. A new test macros has also...
[u/mrichter/AliRoot.git] / ITS / SPDclusterTestBari.C
1 void 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
480 gROOT->Reset();
481 gStyle->SetFillColor(0);
482 gStyle->SetStatW(0.37);
483 gStyle->SetStatH(0.22);
484
485 //----------------------------------------------------- page 1
486 gStyle->SetOptLogy(0);
487 gStyle->SetOptStat(1100);
488
489   TCanvas *c1 = new TCanvas("c1","hits, digits, clusters",200,10,700,780);
490 c1->SetFillColor(0);
491 c1->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 //
609 gStyle->SetOptStat(0);
610 c4->SetFillColor(0);
611 c4->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
675 void 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