ITS/UPGRADE/compression/PatternAnalysis/complete1.C

   1 /*
   2 This macro creates for each pattern four 2D histograms: distance between MChit and COG vs impact angle. This information are taken
   3 from the output of compClusHits.C. Information about the pattern are taken from the database generated by BuildClTop.C.
   4 The macro fits these histograms with gaussians along y axis, in order to study the dependence of dZ and dX wrt
   5 the imact angles. Two 1D histograms with dZ and dX are generated as a projection on to the y axis and are used to calculate mean value-sigma
   6 for each pattern.
   7 The output is a PDF with 4 2D histograms (mean-value of dZ(dX) vs alpha, sigma of dZ(dX) vs alpha,mean-value of dZ(dX) vs beta,
   8 sigma of dZ(dX) vs beta) and 2 histograms, dZ and dX, for each pattern, a drawing of the pattern and the information from the fit and the DB and
   9 information from DB and fits.
  10 You can set the range of patterns you want to study by setting LowerLimit/UpperLimit. If the range coincide with the the whole ragne of patterns,
  11 the DB is updated with the information from these fits.
  12 It is possible to loop over few clusters instead of all the cluster setting ntotCl: usefull to check the macro. (-1 means all the clusters).
  13 It is possilble toprint a PDF with the 2D histograms, pattern drawing and info: setting number2d histos you decide the number of pattern for which
  14 you want 2D histograms.
  15 */
  16
  17
  18 #if !defined(__CINT__) || defined(__MAKECINT__)
  19 #include "TObjArray.h"
  20 #include "TFile.h"
  21 #include "TTree.h"
  22 #include "TH1.h"
  23 #include "TH2.h"
  24 #include "TF1.h"
  25 #include "TArrayI.h"
  26 #include "TArrayF.h"
  27 #include "TCanvas.h"
  28 #include "TLegend.h"
  29 #include "TPad.h"
  30 #include "TLatex.h"
  31 #include "TBits.h"
  32 #include "TStopwatch.h"
  33 #include "TLine.h"
  34 #include "TMarker.h"
  35 #include "../ITS/UPGRADE/AliITSUClusterPix.h"
  36 #include "../ITS/UPGRADE/AliITSURecoLayer.h"
  37 #include "../ITS/UPGRADE/AliITSURecoDet.h"
  38 #include "../ITS/UPGRADE/AliITSUHit.h"
  39 #include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
  40 #include "AliITSsegmentation.h"
  41 #include "AliGeomManager.h"
  42 #include "TROOT.h"
  43 #include "TStyle.h"
  44
  45 #endif
  46
  47 TObjArray* pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
  48 TVectorF* pattFR=0; //frequency of the pattern in the database
  49 TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
  50 TVectorF* zCentrePix=0;
  51 TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
  52 TVectorF* zCentreShift=0;
  53 TVectorF* NPix=0; //Number of fired pixels
  54 TVectorF* NRow=0; //Number of rows of the pattern
  55 TVectorF* NCol=0; //Number of columns of the pattern
  56
  57 TArrayF DeltaZmean; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
  58 TArrayF DeltaZmeanErr;
  59 TArrayF DeltaXmean; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
  60 TArrayF DeltaXmeanErr;
  61 TArrayF DeltaZsigma; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
  62 TArrayF DeltaZsigmaErr;
  63 TArrayF DeltaXsigma; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
  64 TArrayF DeltaXsigmaErr;
  65 TArrayF Chi2x; //Chi2 of the gaussian fit over x
  66 TArrayF Chi2z; //Chi2 of the aussian fut over z
  67 TArrayI NDFx;
  68 TArrayI NDFz;
  69
  70 TObjArray histoArr;
  71
  72 TStopwatch timer;
  73 TStopwatch timer1;
  74
  75 Float_t pitchX; //pitch of the pixel in x direction
  76 Float_t pitchZ; //pitch of the pixel in z direction
  77
  78 enum{kDXAlpha=0, kDZAlpha=1, kDXBeta=2, kDZBeta=3};//to select the kind of TH2 we want
  79
  80 void LoadDB(const char* fname);
  81 void CreateHistos(TObjArray* harr);
  82 TH2* GetHisto(Int_t pattID, Int_t kind, TObjArray* harr);
  83 void FitHistos(Int_t pattID, Float_t Freq, Int_t nPix, Int_t nRow, Int_t nCol,
  84  TObjArray* harr, TCanvas* canvas, TCanvas* canvasbis);
  85 void UpDateDB(const char* outDBFile);
  86
  87 typedef struct {
  88   Int_t evID;
  89   Int_t volID;
  90   Int_t lrID;
  91   Int_t clID;
  92   Int_t nPix;
  93   Int_t nX;
  94   Int_t nZ;
  95   Int_t q;
  96   Float_t pt;
  97   Float_t eta;
  98   Float_t phi;
  99   Float_t xyz[3];
 100   Float_t dX;
 101   Float_t dY;
 102   Float_t dZ;
 103   Bool_t split;
 104   Bool_t prim;
 105   Int_t  pdg;
 106   Int_t  ntr;
 107   Float_t alpha; // alpha is the angle in y-radius plane in local frame
 108   Float_t beta;  // beta is the angle in xz plane, taken from z axis, growing counterclockwise
 109   Int_t nRowPatt;
 110   Int_t nColPatt;
 111   Int_t pattID;
 112   Float_t freq;
 113   Float_t xCen;
 114   Float_t zCen;
 115   Float_t zShift;
 116   Float_t xShift;
 117 } clSumm;
 118
 119 static clSumm Summ;
 120
 121 Int_t nPatterns=0;
 122 Int_t LowerLimit=0;//set the first cluster of the range of interest, INCLUDED (to update DB LowerLimit=0)
 123 Int_t UpperLimit=4000; // EXCLUDED; -1 means number of Patterns. In this case the database is updated with the information about DeltaxMean, DeltaXerr, etc.
 124 Int_t nStudied=0;
 125 Int_t ntotCl=-1;//Set -1 to loop over all the clusters, otherwise set the number of clusters
 126 Int_t number2dhisto=10;//Set the number pattern from 0 to number2dhisto-1 for which you want to print the TH2 histos (dZ or dX vs alpha or beta)
 127
 128 void complete1(){
 129
 130   //importing data
 131
 132   LoadDB("clusterTopology.root");
 133
 134   AliGeomManager::LoadGeometry("geometry.root");
 135   AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
 136   AliITSUClusterPix::SetGeom(gm);
 137   const AliITSsegmentation* segm = gm->GetSegmentation(0);
 138   pitchX = segm->Dpx(0)*10000; //pitch in X in micron
 139   printf("pitchX: %f\n",pitchX);
 140   pitchZ = segm->Dpz(0)*10000; //pitch in Z in micron
 141   printf("pitchX: %f\n",pitchX);
 142   delete gm;
 143
 144   nPatterns = (pattDB->GetEntriesFast());
 145
 146   if(UpperLimit==-1){
 147     UpperLimit=nPatterns;
 148   }
 149
 150   nStudied=UpperLimit-LowerLimit+1;
 151
 152   printf("\n\n nPatterns %d\n\n",nPatterns);
 153
 154   DeltaZmean.Set(nPatterns+100);
 155   DeltaZmeanErr.Set(nPatterns+100);
 156   DeltaXmean.Set(nPatterns+100);
 157   DeltaXmeanErr.Set(nPatterns+100);
 158   DeltaZsigma.Set(nPatterns+100);
 159   DeltaZsigmaErr.Set(nPatterns+100);
 160   DeltaXsigma.Set(nPatterns+100);
 161   DeltaXsigmaErr.Set(nPatterns+100);
 162   Chi2x.Set(nPatterns+100);
 163   Chi2z.Set(nPatterns+100);
 164   NDFx.Set(nPatterns+100);
 165   NDFz.Set(nPatterns+100);
 166
 167   printf("Loading input tree... ");
 168
 169   TFile *input = new TFile("clInfo.root","read");
 170   TTree *clitsu = (TTree*) input->Get("clitsu");
 171
 172   if(ntotCl==-1){
 173     ntotCl = (Int_t) clitsu->GetEntries();
 174   }
 175   TBranch *brdZ = clitsu->GetBranch("dZ");
 176   brdZ->SetAddress(&Summ.dZ);
 177   TBranch *brdX = clitsu->GetBranch("dX");
 178   brdX->SetAddress(&Summ.dX);
 179   TBranch *brpattID = clitsu->GetBranch("pattID");
 180   brpattID->SetAddress(&Summ.pattID);
 181   TBranch *brfreq = clitsu->GetBranch("freq");
 182   brfreq->SetAddress(&Summ.freq);
 183   TBranch *bralpha = clitsu->GetBranch("alpha");
 184   bralpha->SetAddress(&Summ.alpha);
 185   TBranch *brbeta = clitsu->GetBranch("beta");
 186   brbeta->SetAddress(&Summ.beta);
 187   printf("done!!\n\n");
 188
 189   CreateHistos(&histoArr);
 190   printf("Starting to process clusters\n\n");
 191   timer.Start();
 192   for(Int_t iCl=0; iCl<ntotCl; iCl++){
 193     clitsu->GetEntry(iCl);
 194     Int_t ID = Summ.pattID;
 195     printf("Processing cluster %d... ",iCl);
 196     if(ID<LowerLimit || ID>=UpperLimit) {
 197       printf("skipped!\n\n");
 198       continue;
 199     }
 200     GetHisto(ID,kDXAlpha,&histoArr)->Fill(Summ.alpha,Summ.dX);
 201     GetHisto(ID,kDZAlpha,&histoArr)->Fill(Summ.alpha,Summ.dZ);
 202     GetHisto(ID,kDXBeta,&histoArr)->Fill(Summ.beta,Summ.dX);
 203     GetHisto(ID,kDZBeta,&histoArr)->Fill(Summ.beta,Summ.dZ);
 204     printf("done!\n\n");
 205   }
 206   printf("All clusters processed!\n");
 207   timer.Stop();
 208   timer.Print();
 209
 210   TCanvas* cnv1 = new TCanvas("cnv1","cnv1",900,700);
 211   TCanvas* cnv2 = new TCanvas("cnv2","cnv2",900,700);
 212
 213   printf("Printing PDF...\n\n");
 214   timer1.Start();
 215
 216   cnv1->Print("angles.pdf[");
 217   cnv2->Print("angles2D.pdf[");
 218
 219   for(Int_t i=LowerLimit; i<UpperLimit; i++){
 220
 221     FitHistos(i,(*pattFR)[i],(*NPix)[i],(*NRow)[i],(*NCol)[i],&histoArr, cnv1, cnv2);
 222   }
 223
 224   cnv1->Print("angles.pdf]");
 225   cnv2->Print("angles2D.pdf]");
 226   timer1.Stop();
 227   printf("\n\ndone!!\n\n");
 228   timer1.Print();
 229
 230   delete cnv1;
 231   delete cnv2;
 232
 233   if(LowerLimit==0 && UpperLimit==nPatterns){
 234     UpDateDB("clusterTopology2.root");
 235   }
 236
 237 }
 238
 239 void LoadDB(const char* fname){ //load the data base built by BuildClTopDB.C
 240
 241   // load database
 242   TFile* fl = TFile::Open(fname);
 243   if(!fl){printf("Could not find %s",fname); exit(1);}
 244   pattDB = (TObjArray*)fl->Get("TopDB");
 245   pattFR = (TVectorF*)fl->Get("TopFreq");
 246   xCentrePix =(TVectorF*)fl->Get("xCOG");
 247   zCentrePix =(TVectorF*)fl->Get("zCOG");
 248   xCentreShift =(TVectorF*)fl->Get("xShift");
 249   zCentreShift =(TVectorF*)fl->Get("zShift");
 250   NPix =(TVectorF*)fl->Get("NPix");
 251   NCol =(TVectorF*)fl->Get("NCol");
 252   NRow =(TVectorF*)fl->Get("NRow");
 253 }
 254
 255 void CreateHistos(TObjArray* harr){ //creates four 2D histograms (dX or dZ vs alpha or beta) for each pattern included in the study range
 256
 257   printf("Creating histograms... ");
 258
 259   //for each pattern we crate 4 2D histograms corresponding the cathegories within the enumeration
 260   if(!harr) harr=&histoArr;
 261
 262   for(Int_t i=0; i<nStudied; i++){
 263
 264     TH2* h0 = new TH2F(Form("hXalpha%d", i),
 265       Form("#DeltaX vs #alpha for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
 266     h0->SetDirectory(0);
 267     harr->AddAtAndExpand(h0, i*4);//This histogram is the first of the 4-histos series
 268
 269     TH2* h1 = new TH2F(Form("hZalpha%d", i),
 270       Form("#DeltaZ vs #alpha for pattID %d", i),10, -0.1, 1.1* TMath::Pi()/2,100,-30,30);
 271     h1->SetDirectory(0);
 272     harr->AddAtAndExpand(h1, i*4+1);//This histogram is the second of the 4-histos series and so on...
 273
 274     TH2* h2 = new TH2F(Form("hXbeta%d", i),
 275       Form("#DeltaX vs #beta for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
 276     h2->SetDirectory(0);
 277     harr->AddAtAndExpand(h2, i*4+2);
 278
 279     TH2* h3 = new TH2F(Form("hZbeta%d", i),
 280       Form("#DeltaZ vs #beta for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
 281     h3->SetDirectory(0);
 282     harr->AddAtAndExpand(h3, i*4+3);
 283
 284   }
 285
 286   printf("done!!\n\n");
 287 }
 288
 289 TH2* GetHisto(Int_t pattID, Int_t kind, TObjArray* harr){
 290
 291   if(!harr) harr=&histoArr ;
 292   TH2*  h=0;
 293   h=(TH2*)harr->At((pattID-LowerLimit)*4+kind);
 294   if (!h) {
 295     printf("Unknown histo id=%d\n",pattID);
 296     exit(1);
 297   }
 298   return h;
 299
 300 }
 301
 302 void FitHistos(Int_t pattID, Float_t Freq, Int_t nPix, Int_t nRow,
 303  Int_t nCol, TObjArray* harr, TCanvas* canvas, TCanvas* canvasbis){ //Fit the 2D histograms with gaussian along the Y and creates two 1D histograms,
 304                                                                     // with dX and dZ, bay a projection to the y axis and the fit with gaussian
 305
 306   static TF1* gs = new TF1("gs","gaus",-50,50);
 307
 308   if(!harr) harr=&histoArr;
 309
 310   //Creating the pads for the output
 311
 312   canvas->Clear();
 313   canvas->cd();
 314
 315   Float_t width = (1 - canvas->GetRightMargin() - canvas->GetLeftMargin() )/3;
 316
 317   TPad* titlepad = new TPad("titlepad","",canvas->GetLeftMargin(),0.825,1-canvas->GetRightMargin(),0.95);
 318   titlepad->SetFillColor(kYellow);
 319   titlepad->Draw();
 320   TPad* pad1 = new TPad("pad1","",canvas->GetLeftMargin(),(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),canvas->GetLeftMargin()+width,0.8);
 321   pad1->Draw();
 322   TPad* pad2 = new TPad("pad2","",canvas->GetLeftMargin()+width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),canvas->GetLeftMargin()+2*width,0.8);
 323   pad2->Draw();
 324   TPad* pad3 = new TPad("pad3","",canvas->GetLeftMargin()+2*width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),1-canvas->GetRightMargin(),0.8);
 325   pad3->Draw();
 326   TPad* pad4 = new TPad("pad3","",canvas->GetLeftMargin(),canvas->GetBottomMargin(),canvas->GetLeftMargin()+width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
 327   pad4->Draw();
 328   TPad* pad5 = new TPad("pad5","",canvas->GetLeftMargin()+width,canvas->GetBottomMargin(),canvas->GetLeftMargin()+2*width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
 329   pad5->Draw();
 330   TPad* pad6 = new TPad("pad6","",canvas->GetLeftMargin()+2*width,canvas->GetBottomMargin(),1-canvas->GetRightMargin(),(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
 331   pad6->Draw();
 332   TPad* padPattern = new TPad("padPattern","",0.75,0.825,0.9,0.95);
 333   padPattern->Draw();
 334
 335   //...................................................................................................
 336
 337   // dX vs alpha
 338
 339   TObjArray harrFXA;
 340   TH2* hXA = (TH2*)harr->At((pattID-LowerLimit)*4);
 341   hXA->GetYaxis()->SetTitle("#DeltaX, #mum");
 342   hXA->GetXaxis()->SetTitle("#alpha");
 343   hXA->SetStats(0);
 344   hXA->FitSlicesY(gs,0,-1,0,"qnr",&harrFXA);
 345   TH1* hXA_1 = (TH1*)harrFXA.At(1);
 346   hXA_1->SetStats(0);
 347   hXA_1->SetDirectory(0);
 348   hXA_1->GetYaxis()->SetTitle("#DeltaX_{#mu}(#DeltaZ_{#mu}), #mum");
 349   TH1* hXA_2 = (TH1*)harrFXA.At(2);
 350   hXA_2->SetStats(0);
 351   hXA_2->SetDirectory(0);
 352   hXA_2->GetYaxis()->SetTitle("#DeltaX_{#sigma}(#DeltaZ_{#sigma}), #mum");
 353
 354   //............................................................................................................
 355
 356   // dX
 357
 358   Int_t fitStatusX=0;
 359
 360   TH1* hdx = hXA->ProjectionY("hdx");
 361   hdx->GetXaxis()->SetTitle("#DeltaX, #mum");
 362   hdx->SetTitle(Form("#DeltaX for pattID %d",pattID));
 363   hdx->SetStats(0);
 364   pad3->cd();
 365   gs->SetParameters(hdx->GetMaximum(),hdx->GetMean(),hdx->GetRMS());
 366   if((hdx->GetEntries())<100) fitStatusX = hdx->Fit("gs","ql");
 367   else fitStatusX = hdx->Fit("gs","q");
 368
 369   if(fitStatusX==0){
 370     DeltaXmean[pattID]= gs->GetParameter(1);
 371     DeltaXmeanErr[pattID]= gs->GetParError(1);
 372     DeltaXsigma[pattID]= gs->GetParameter(2);
 373     DeltaXsigmaErr[pattID]= gs->GetParError(2);
 374     Chi2x[pattID] = gs->GetChisquare();
 375     Int_t varNDFx = gs->GetNDF();
 376     if(varNDFx>=0)
 377     NDFx[pattID] = varNDFx;
 378     else{
 379       NDFx[pattID]=-1;
 380     }
 381   }
 382   else{
 383     DeltaXmean[pattID]=0;
 384     DeltaXmeanErr[pattID]=0;
 385     DeltaXsigma[pattID]=0;
 386     DeltaXsigmaErr[pattID]=0;
 387     Chi2x[pattID] = -1;
 388   }
 389
 390   pad3->cd();
 391   hdx->Draw();
 392
 393
 394   //.............................................................................................................
 395
 396
 397   // dZ vs alpha
 398
 399   TObjArray harrFZA;
 400   TH2* hZA = (TH2*)harr->At((pattID-LowerLimit)*4+1);
 401   hZA->GetYaxis()->SetTitle("#DeltaZ, #mum");
 402   hZA->GetXaxis()->SetTitle("#alpha");
 403   hZA->SetStats(0);
 404   hZA->FitSlicesY(gs,0,-1,0,"qnr",&harrFZA);
 405   TH1* hZA_1 = (TH1*)harrFZA.At(1);
 406   hZA_1->SetMarkerColor(8);
 407   hZA_1->SetLineColor(8);
 408   hZA_1->SetStats(0);
 409   TH1* hZA_2 = (TH1*)harrFZA.At(2);
 410   hZA_2->SetMarkerColor(8);
 411   hZA_2->SetLineColor(8);
 412   hZA_2->SetStats(0);
 413
 414   Double_t maxA_1 = TMath::Max(hXA_1->GetMaximum(),hZA_1->GetMaximum());
 415   Double_t minA_1 = TMath::Min(hXA_1->GetMinimum(),hZA_1->GetMinimum());
 416   Double_t addtorangeA_1 = (maxA_1-minA_1)*0.1;
 417   Double_t rangeA_1 = addtorangeA_1*10;
 418
 419   Double_t maxA_2 = TMath::Max(hXA_2->GetMaximum(),hZA_2->GetMaximum());
 420   Double_t minA_2 = TMath::Min(hXA_2->GetMinimum(),hZA_2->GetMinimum());
 421   Double_t addtorangeA_2 = (maxA_2-minA_2)*0.1;
 422   Double_t rangeA_2 = addtorangeA_2*10;
 423
 424   hXA_1->SetMaximum(maxA_1+addtorangeA_1);
 425   hXA_1->SetMinimum(minA_1-addtorangeA_1);
 426
 427   hXA_2->SetMaximum(maxA_2+addtorangeA_2);
 428   hXA_2->SetMinimum(minA_2-addtorangeA_2);
 429
 430   pad1->cd();
 431   hXA_1->Draw();
 432   hZA_1->Draw("same");
 433   TLegend* legA_1 = new TLegend(0.8*1.1* TMath::Pi()/2, maxA_1-0.2*rangeA_1, TMath::Pi()/2, maxA_1,"","");
 434   legA_1->AddEntry(hXA_1, "   #DeltaX_{#mu}", "l");
 435   legA_1->AddEntry(hZA_1, "   #DeltaZ_{#mu}", "l");
 436   legA_1->SetFillColor(0);
 437   legA_1->SetBorderSize(0);
 438   legA_1->Draw();
 439   pad2->cd();
 440   hXA_2->Draw();
 441   hZA_2->Draw("same");
 442   TLegend* legA_2 = new TLegend(0.8*1.1* TMath::Pi()/2, maxA_2-0.2*rangeA_2, TMath::Pi()/2, maxA_2,"","");
 443   legA_2->AddEntry(hXA_2, "   #DeltaX_{#sigma}", "l");
 444   legA_2->AddEntry((TObject*)0, "", "");
 445   legA_2->AddEntry(hZA_2, "   #DeltaZ_{#sigma}", "l");
 446   legA_2->SetFillColor(0);
 447   legA_2->SetBorderSize(0);
 448   legA_2->Draw();
 449
 450   //...........................................................................................
 451
 452
 453   // dZ
 454
 455   Int_t fitStatusZ=0;
 456
 457   TH1* hdz = hZA->ProjectionY("hdz");
 458   hdz->GetXaxis()->SetTitle("#DeltaZ, #mum");
 459   hdz->SetTitle(Form("#DeltaZ for pattID %d",pattID));
 460   hdz->SetStats(0);
 461   pad6->cd();
 462   gs->SetParameters(hdz->GetMaximum(),hdz->GetMean(),hdz->GetRMS());
 463   if((hdz->GetEntries())<100) fitStatusZ = hdz->Fit("gs","ql");
 464   else fitStatusZ = hdz->Fit("gs","q");
 465
 466   if(fitStatusZ==0){
 467     DeltaZmean[pattID]= gs->GetParameter(1);
 468     DeltaZmeanErr[pattID]= gs->GetParError(1);
 469     DeltaZsigma[pattID]= gs->GetParameter(2);
 470     DeltaZsigmaErr[pattID]= gs->GetParError(2);
 471     Chi2z[pattID] = gs->GetChisquare();
 472     Int_t varNDFz = gs->GetNDF();
 473     if(varNDFz>=0)
 474     NDFz[pattID] = varNDFz;
 475     else{
 476     NDFz[pattID]=-1;
 477     }
 478   }
 479   else{
 480     DeltaZmean[pattID]=0;
 481     DeltaZmeanErr[pattID]=0;
 482     DeltaZsigma[pattID]=0;
 483     DeltaZsigmaErr[pattID]=0;
 484     Chi2z[pattID] = -1;
 485   }
 486
 487   pad6->cd();
 488   hdz->Draw();
 489
 490   //.............................................................................................................
 491
 492   // dX vs beta
 493
 494   TObjArray harrFXB;
 495   TH2* hXB = (TH2*)harr->At((pattID-LowerLimit)*4+2);
 496   hXB->GetYaxis()->SetTitle("#DeltaX, #mum");
 497   hXB->GetXaxis()->SetTitle("#beta");
 498   hXB->SetStats(0);
 499   hXB->FitSlicesY(gs,0,-1,0,"qnr",&harrFXB);
 500   TH1* hXB_1 = (TH1*)harrFXB.At(1);
 501   hXB_1->SetDirectory(0);
 502   hXB_1->SetStats(0);
 503   hXB_1->GetYaxis()->SetTitle("#DeltaX_{#mu}(#DeltaZ_{#mu}), #mum");
 504   TH1* hXB_2 = (TH1*)harrFXB.At(2);
 505   hXB_2->SetDirectory(0);
 506   hXB_2->SetStats(0);
 507   hXB_2->GetYaxis()->SetTitle("#DeltaX_{#sigma}(#DeltaZ_{#sigma}), #mum");
 508
 509   //............................................................................................
 510
 511   //dZ vs beta
 512
 513   TObjArray harrFZB;
 514   TH2* hZB = (TH2*)harr->At((pattID-LowerLimit)*4+3);
 515   hZB->GetYaxis()->SetTitle("#DeltaZ, #mum");
 516   hZB->GetXaxis()->SetTitle("#alpha");
 517   hZB->SetStats(0);
 518   hZB->FitSlicesY(gs,0,-1,0,"qnr",&harrFZB);
 519   TH1* hZB_1 = (TH1*)harrFZB.At(1);
 520   hZB_1->SetMarkerColor(8);
 521   hZB_1->SetLineColor(8);
 522   hZB_1->SetStats(0);
 523   TH1* hZB_2 = (TH1*)harrFZB.At(2);
 524   hZB_2->SetMarkerColor(8);
 525   hZB_2->SetLineColor(8);
 526   hZB_2->SetStats(0);
 527
 528   Double_t maxB_1 = TMath::Max(hXB_1->GetMaximum(),hZB_1->GetMaximum());
 529   Double_t minB_1 = TMath::Min(hXB_1->GetMinimum(),hZB_1->GetMinimum());
 530   Double_t addtorangeB_1 = (maxB_1-minB_1)*0.1;
 531   Double_t rangeB_1 = 10*addtorangeB_1;
 532
 533   Double_t maxB_2 = TMath::Max(hXB_2->GetMaximum(),hZB_2->GetMaximum());
 534   Double_t minB_2 = TMath::Min(hXB_2->GetMinimum(),hZB_2->GetMinimum());
 535   Double_t addtorangeB_2 = (maxB_2-minB_2)*0.1;
 536   Double_t rangeB_2 = 10*addtorangeB_2;
 537
 538
 539   hXB_1->SetMaximum(maxB_1+addtorangeB_1);
 540   hXB_1->SetMinimum(minB_1-addtorangeB_1);
 541
 542   hXB_2->SetMaximum(maxB_2+addtorangeB_2);
 543   hXB_2->SetMinimum(minB_2-addtorangeB_2);
 544
 545   pad4->cd();
 546   hXB_1->Draw();
 547   hZB_1->Draw("same");
 548   TLegend* legB_1 = new TLegend(0.8*1.1* TMath::Pi()/2, maxB_1-0.2*rangeB_1, TMath::Pi()/2, maxB_1,"","");
 549   legB_1->AddEntry(hXA_1, "   #DeltaX_{#mu}", "l");
 550   legB_1->AddEntry((TObject*)0, "", "");
 551   legB_1->AddEntry(hZA_1, "   #DeltaZ_{#mu}", "l");
 552   legB_1->SetFillColor(0);
 553   legB_1->SetBorderSize(0);
 554   legB_1->Draw();
 555   pad5->cd();
 556   hXB_2->Draw();
 557   hZB_2->Draw("same");
 558   TLegend* legB_2 = new TLegend(0.8*1.1* TMath::Pi()/2, maxB_2-0.2*rangeB_2, TMath::Pi()/2, maxB_2,"","");
 559   legB_2->AddEntry(hXB_2, "   #DeltaX_{#sigma}", "l");
 560   legB_2->AddEntry((TObject*)0, "", "");
 561   legB_2->AddEntry(hZB_2, "   #DeltaZ_{#sigma}", "l");
 562   legB_2->SetFillColor(0);
 563   legB_2->SetBorderSize(0);
 564   legB_2->Draw();
 565
 566   //cnv->Print("patt0.gif");
 567
 568   //..................................................................................................................
 569
 570   // Drawing the title
 571
 572   const char* text = Form("pattID: %d   freq: %f   nPix: %d   nRow: %d   nCol: %d  xCOG (fraction): %0.2f + (%0.2f)   zCOG (fraction): %0.2f + (%0.2f)",
 573     pattID,Freq,nPix,nRow,nCol, (*xCentrePix)[pattID],(*xCentreShift)[pattID],(*zCentrePix)[pattID],
 574     (*zCentreShift)[pattID]);
 575   const char* text1 = Form("#DeltaX_{#mu}: (%.2f #pm %.2f) #mum     #DeltaX_{#sigma}: (%.2f #pm %.2f) #mum     #DeltaX_{MC-centre}: %.2f #mum     #chi^{2}/NDF: %.2f / %d",
 576     DeltaXmean[pattID],DeltaXmeanErr[pattID],DeltaXsigma[pattID],DeltaXsigmaErr[pattID], DeltaXmean[pattID]+(*xCentreShift)[pattID]*pitchX, Chi2x[pattID], NDFx[pattID]);
 577   const char* text2 = Form("#DeltaZ_{#mu}: (%.2f #pm %.2f) #mum     #DeltaZ_{#sigma}: (%.2f #pm %.2f) #mum     #DeltaZ_{MC-centre}: %.2f #mum     #chi^{2}/NDF: %.2f / %d",
 578     DeltaZmean[pattID],DeltaZmeanErr[pattID],DeltaZsigma[pattID],DeltaZsigmaErr[pattID], DeltaZmean[pattID]+(*zCentreShift)[pattID]*pitchZ, Chi2z[pattID], NDFz[pattID]);
 579   TLatex title;
 580   title.SetTextSize(0.1);
 581   title.SetTextAlign(12);
 582   titlepad->cd();
 583   title.DrawLatex(0.085,0.75,text);
 584   title.DrawLatex(0.085,0.5,text1);
 585   title.DrawLatex(0.085,0.25,text2);
 586
 587
 588   //Drawing the pattern
 589
 590   TH2* patternIm = new TH2F("patternIm","", nCol,-0.5,nCol-0.5,nRow,-0.5,nRow-0.5);
 591   for (Int_t ir=0; ir<nRow; ir++) {
 592     for (Int_t ic=0; ic<nCol; ic++) {
 593       if(((TBits*)(*pattDB)[pattID])->TestBitNumber(ir*nCol+ic)){
 594         patternIm->Fill(ic,ir);
 595       }
 596     }
 597   }
 598
 599   patternIm->SetStats(0);
 600   padPattern->cd();
 601   patternIm->Draw("Acol");
 602
 603   //Drawing vertical lines on the pattern
 604
 605   TObjArray vertlinesarray;
 606   for(int i=2; i<=nCol; i++){
 607     TLine* vline = new TLine(patternIm->GetXaxis()->GetBinLowEdge(i),patternIm->GetYaxis()->GetBinLowEdge(1),
 608       patternIm->GetXaxis()->GetBinLowEdge(i),patternIm->GetYaxis()->GetBinLowEdge(nRow+1));
 609     vertlinesarray.AddAtAndExpand(vline,i);
 610   }
 611
 612   for(int i=2; i<=nCol; i++){
 613     TLine* v1 = (TLine*)vertlinesarray.At(i);
 614     padPattern->cd();
 615     v1->Draw();
 616   }
 617
 618   //Drawing horizontal lines on the pattern
 619
 620   TObjArray horlinesarray;
 621   for(int i=2; i<=nRow; i++){
 622     TLine* hline = new TLine(patternIm->GetXaxis()->GetBinLowEdge(1),patternIm->GetYaxis()->GetBinLowEdge(i),
 623       patternIm->GetXaxis()->GetBinLowEdge(nCol+1),patternIm->GetYaxis()->GetBinLowEdge(i));
 624     horlinesarray.AddAtAndExpand(hline,i);
 625   }
 626
 627   for(int i=2; i<=nRow; i++){
 628     TLine* hor1 = (TLine*)horlinesarray.At(i);
 629     padPattern->cd();
 630     hor1->Draw();
 631   }
 632
 633   //Drawing MChit and COG
 634
 635   TMarker* COG = new TMarker((*zCentrePix)[pattID]+(*zCentreShift)[pattID]-0.5,
 636     (*xCentrePix)[pattID]+(*xCentreShift)[pattID]-0.5,20);
 637   COG->Draw();
 638
 639   TMarker* MChit = new TMarker((*zCentrePix)[pattID]+(*zCentreShift)[pattID]+(DeltaZmean[pattID])/pitchZ-0.5, //20 is the pixel pitch
 640     (*xCentrePix)[pattID]+(*xCentreShift)[pattID]+(DeltaXmean[pattID])/pitchX-0.5,29);
 641   MChit->SetMarkerColor(kWhite);
 642   MChit->Draw();
 643
 644   canvas->Print("angles.pdf");
 645   canvas->Clear();
 646
 647   if(pattID<number2dhisto){
 648
 649     canvasbis->cd();
 650     TPad* titlepadbis = new TPad("titlepadbis","",canvas->GetLeftMargin(),0.825,1-canvas->GetRightMargin(),0.95);
 651     titlepadbis->SetFillColor(kYellow);
 652     titlepadbis->Draw();
 653     TPad* pad1bis = new TPad("pad1bis","",canvasbis->GetLeftMargin(),(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin(),0.5,0.8);
 654     pad1bis->Draw();
 655     TPad* pad2bis = new TPad("pad2bis","",0.5,(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin(),1-canvasbis->GetRightMargin(),0.8);
 656     pad2bis->Draw();
 657     TPad* pad3bis = new TPad("pad3bis","",canvasbis->GetLeftMargin(),canvasbis->GetBottomMargin(),0.5,(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin());
 658     pad3bis->Draw();
 659     TPad* pad4bis = new TPad("pad4bis","",0.5,canvasbis->GetBottomMargin(),1-canvasbis->GetRightMargin(),(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin());
 660     pad4bis->Draw();
 661     TPad* padPatternbis = new TPad("padPatternbis","",0.75,0.825,0.9,0.95);
 662     padPatternbis->Draw();
 663
 664     titlepadbis->cd();
 665     title.DrawLatex(0.085,0.75,text);
 666     title.DrawLatex(0.085,0.5,text1);
 667     title.DrawLatex(0.085,0.25,text2);
 668     pad1bis->cd();
 669     hXA->Draw("colz");
 670     pad2bis->cd();
 671     hXB->Draw("colz");
 672     pad3bis->cd();
 673     hZA->Draw("colz");
 674     pad4bis->cd();
 675     hZB->Draw("colz");
 676     padPatternbis->cd();
 677     patternIm->Draw("Acol");
 678     for(int i=2; i<=nCol; i++){
 679       TLine* v1 = (TLine*)vertlinesarray.At(i);
 680       padPatternbis->cd();
 681       v1->Draw();
 682     }
 683     for(int i=2; i<=nRow; i++){
 684       TLine* hor1 = (TLine*)horlinesarray.At(i);
 685       padPatternbis->cd();
 686       hor1->Draw();
 687     }
 688     COG->Draw();
 689     MChit->Draw();
 690     canvasbis->Print("angles2D.pdf");
 691     canvasbis->Clear();
 692   }
 693
 694   delete patternIm;
 695   delete COG;
 696   delete MChit;
 697
 698   horlinesarray.Delete();
 699   vertlinesarray.Delete();
 700
 701 }
 702
 703 void UpDateDB(const char* outDBFile){
 704
 705   printf("\n\n\nStoring data in the DataBase\n\n\n");
 706
 707   TString outDBFileName = outDBFile;
 708   if (outDBFileName.IsNull()) outDBFileName = "clusterTopology2.root";
 709
 710   TVectorF arrDeltaZmean(nPatterns);
 711   TVectorF arrDeltaZmeanErr(nPatterns);
 712   TVectorF arrDeltaXmean(nPatterns);
 713   TVectorF arrDeltaXmeanErr(nPatterns);
 714   TVectorF arrDeltaZsigma(nPatterns);
 715   TVectorF arrDeltaZsigmaErr(nPatterns);
 716   TVectorF arrDeltaXsigma(nPatterns);
 717   TVectorF arrDeltaXsigmaErr(nPatterns);
 718   TVectorF arrChi2x(nPatterns);
 719   TVectorF arrNDFx(nPatterns);
 720   TVectorF arrChi2z(nPatterns);
 721   TVectorF arrNDFz(nPatterns);
 722
 723   for(Int_t ID=0; ID<nPatterns; ID++){
 724
 725     printf("Processing pattern %d... ", ID);
 726
 727     arrDeltaZmean[ID]=DeltaZmean[ID];
 728     arrDeltaZmeanErr[ID]=DeltaZmeanErr[ID];
 729     arrDeltaXmean[ID]=DeltaXmean[ID];
 730     arrDeltaXmeanErr[ID]=DeltaXmeanErr[ID];
 731     arrDeltaZsigma[ID]=DeltaZsigma[ID];
 732     arrDeltaZsigmaErr[ID]=DeltaZsigmaErr[ID];
 733     arrDeltaXsigma[ID]=DeltaXsigma[ID];
 734     arrDeltaXsigmaErr[ID]=DeltaXsigmaErr[ID];
 735     arrChi2x[ID]=Chi2x[ID];
 736     arrNDFx[ID]=NDFx[ID];
 737     arrChi2z[ID]=Chi2z[ID];
 738     arrNDFz[ID]=NDFz[ID];
 739
 740     printf("done!!\n\n");
 741   }
 742
 743   TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
 744   flDB->WriteObject(pattDB,"TopDB","kSingleKey");
 745   flDB->WriteObject(NPix,"NPix","kSingleKey");
 746   flDB->WriteObject(NRow,"NRow","kSingleKey");
 747   flDB->WriteObject(NCol,"NCol","kSingleKey");
 748   flDB->WriteObject(pattFR, "TopFreq","kSingleKey");
 749   flDB->WriteObject(xCentrePix,"xCOG","kSingleKey");
 750   flDB->WriteObject(xCentreShift,"xShift","kSingleKey");
 751   flDB->WriteObject(zCentrePix,"zCOG","kSingleKey");
 752   flDB->WriteObject(zCentreShift,"zShift","kSingleKey");
 753   flDB->WriteObject(&arrDeltaZmean,"DeltaZmean","kSingleKey");
 754   flDB->WriteObject(&arrDeltaZmeanErr,"DeltaZmeanErr","kSingleKey");
 755   flDB->WriteObject(&arrDeltaXmean,"DeltaXmean","kSingleKey");
 756   flDB->WriteObject(&arrDeltaXmeanErr,"DeltaXmeanErr","kSingleKey");
 757   flDB->WriteObject(&arrDeltaZsigma,"DeltaZsigma","kSingleKey");
 758   flDB->WriteObject(&arrDeltaZsigmaErr,"DeltaZsigmaErr","kSingleKey");
 759   flDB->WriteObject(&arrDeltaXsigma,"DeltaXsigma","kSingleKey");
 760   flDB->WriteObject(&arrDeltaXsigmaErr,"DeltaXsigmaErr","kSingleKey");
 761   flDB->WriteObject(&arrChi2x,"Chi2x","kSingleKey");
 762   flDB->WriteObject(&arrChi2z,"Chi2z","kSingleKey");
 763   flDB->WriteObject(&arrNDFx,"NDFx","kSingleKey");
 764   flDB->WriteObject(&arrNDFz,"NDFz","kSingleKey");
 765   //
 766   flDB->Close();
 767   delete flDB;
 768
 769   printf("\n\nDB Complete!!\n\n");
 770 }