ITS/SPDclusterTest.C

   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
  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