ITS/UPGRADE/compression/PatternAnalysis/errClus.C

   1 /*
   2 This macro creates for each pattern in the database (generated by BuildClTop.C) two histograms, for MChit-COG difference in x and z direction respectively
   3 This histograms are filled with the data from compClusHits.C. With a gaussian fit we get mean value, sigma, chi2, NDF and update the database of patterns
   4 with this information.
   5 */
   6
   7 #if !defined(__CINT__) || defined(__MAKECINT__)
   8 #include "TObjArray.h"
   9 #include "TFile.h"
  10 #include "TTree.h"
  11 #include "TH1F.h"
  12 #include "TH2F.h"
  13 #include "TArrayI.h"
  14 #include "TArrayF.h"
  15 #include "TCanvas.h"
  16 #include "TBits.h"
  17 #include "TF1.h"
  18 #include "../ITS/UPGRADE/AliITSUClusterPix.h"
  19 #include "../ITS/UPGRADE/AliITSURecoLayer.h"
  20 #include "../ITS/UPGRADE/AliITSURecoDet.h"
  21 #include "../ITS/UPGRADE/AliITSUHit.h"
  22 #include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
  23 #include "AliITSsegmentation.h"
  24 #include "AliGeomManager.h"
  25 #include "TROOT.h"
  26 #include "TStyle.h"
  27 #include "TCanvas.h"
  28 #include "TPaveStats.h"
  29
  30 #endif
  31
  32 TObjArray *pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
  33 TVectorF* pattFR=0; //frequecy of the pattern in the database
  34 TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
  35 TVectorF* zCentrePix=0;
  36 TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
  37 TVectorF* zCentreShift=0;
  38 TVectorF* NPix=0; //Number of fired pixels
  39 TVectorF* NRow=0; //Number of rows of the pattern
  40 TVectorF* NCol=0; //Number of columns of the pattern
  41
  42 TArrayF DeltaZmean; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
  43 TArrayF DeltaZmeanErr;
  44 TArrayF DeltaXmean; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
  45 TArrayF DeltaXmeanErr;
  46 TArrayF DeltaZsigma; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
  47 TArrayF DeltaZsigmaErr;
  48 TArrayF DeltaXsigma; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
  49 TArrayF DeltaXsigmaErr;
  50 TArrayF Chi2x; //Chi2 of the gaussian fit over x
  51 TArrayF Chi2z; //Chi2 of the aussian fut over z
  52 TArrayI NDFx;
  53 TArrayI NDFz;
  54
  55 TObjArray histoArr;
  56
  57 /*Int_t err1counter=0;
  58 Int_t err2counter=0;
  59 Int_t err3counter=0;*/
  60
  61 enum{kXhisto=0,kZhisto=1};
  62
  63 void LoadDB(const char* fname);
  64 void CreateHistos(TObjArray* harr);
  65 TH1* GetHistoX(Int_t pattID, TObjArray* harr);
  66 TH1* GetHistoZ(Int_t pattID, TObjArray* harr);
  67 void DrawReport(TObjArray* harr);
  68 void DrawNP(int pattID, TObjArray* harr);
  69 void UpDateDB(const char* outDBFile);
  70
  71 typedef struct {
  72   Int_t evID;
  73   Int_t volID;
  74   Int_t lrID;
  75   Int_t clID;
  76   Int_t nPix;
  77   Int_t nX;
  78   Int_t nZ;
  79   Int_t q;
  80   Float_t pt;
  81   Float_t eta;
  82   Float_t phi;
  83   Float_t xyz[3];
  84   Float_t dX;
  85   Float_t dY;
  86   Float_t dZ;
  87   Bool_t split;
  88   Bool_t prim;
  89   Int_t  pdg;
  90   Int_t  ntr;
  91   Float_t alpha; // alpha is the angle in y-radius plane in local frame
  92   Float_t beta;  // beta is the angle in xz plane, taken from z axis, growing counterclockwise
  93   Int_t nRowPatt;
  94   Int_t nColPatt;
  95   Int_t pattID;
  96   Float_t freq;
  97   Float_t xCen;
  98   Float_t zCen;
  99   Float_t zShift;
 100   Float_t xShift;
 101 } clSumm;
 102
 103 static clSumm Summ;
 104
 105 Int_t nPatterns=0;
 106
 107 void errClus(){
 108
 109         //importing data
 110
 111   LoadDB("clusterTopology.root");
 112
 113   nPatterns = (pattDB->GetEntriesFast());
 114
 115   printf("\n\n nPatterns %d\n\n",nPatterns);
 116
 117   DeltaZmean.Set(nPatterns+100);
 118   DeltaZmeanErr.Set(nPatterns+100);
 119   DeltaXmean.Set(nPatterns+100);
 120   DeltaXmeanErr.Set(nPatterns+100);
 121   DeltaZsigma.Set(nPatterns+100);
 122   DeltaZsigmaErr.Set(nPatterns+100);
 123   DeltaXsigma.Set(nPatterns+100);
 124   DeltaXsigmaErr.Set(nPatterns+100);
 125   Chi2x.Set(nPatterns+100);
 126   Chi2z.Set(nPatterns+100);
 127   NDFx.Set(nPatterns+100);
 128   NDFz.Set(nPatterns+100);
 129
 130   TFile *input = new TFile("clInfo.root","read");
 131         TTree *clitsu = (TTree*) input->Get("clitsu");
 132
 133         Int_t ntotCl = (Int_t) clitsu->GetEntries();
 134         TBranch *brdZ = clitsu->GetBranch("dZ");
 135         brdZ->SetAddress(&Summ.dZ);
 136   TBranch *brdX = clitsu->GetBranch("dX");
 137   brdX->SetAddress(&Summ.dX);
 138   TBranch *brpattID = clitsu->GetBranch("pattID");
 139   brpattID->SetAddress(&Summ.pattID);
 140
 141   CreateHistos(&histoArr);
 142
 143   for(Int_t iCl=0; iCl<ntotCl; iCl++){
 144     clitsu->GetEntry(iCl);
 145     GetHistoX(Summ.pattID, &histoArr)->Fill(Summ.dX);
 146     GetHistoZ(Summ.pattID, &histoArr)->Fill(Summ.dZ);
 147   }
 148
 149   DrawReport(&histoArr);
 150
 151   //printf("\n\nNDF is 0 : %d times\n\nNDF and fitstati not corresponding: %d times\n\nNDF not corresponding : %d times\n\n", err3counter,err1counter, err2counter);
 152
 153   UpDateDB("clusterTopology2.root");
 154
 155 }
 156
 157 void LoadDB(const char* fname){ //load the data base built by BuildClTopDB.C
 158
 159   // load database
 160   TFile* fl = TFile::Open(fname);
 161   if(!fl){printf("Could not find %s",fname); exit(1);}
 162   pattDB = (TObjArray*)fl->Get("TopDB");
 163   pattFR = (TVectorF*)fl->Get("TopFreq");
 164   xCentrePix =(TVectorF*)fl->Get("xCOG");
 165   zCentrePix =(TVectorF*)fl->Get("zCOG");
 166   xCentreShift =(TVectorF*)fl->Get("xShift");
 167   zCentreShift =(TVectorF*)fl->Get("zShift");
 168   NPix =(TVectorF*)fl->Get("NPix");
 169   NCol =(TVectorF*)fl->Get("NCol");
 170   NRow =(TVectorF*)fl->Get("NRow");
 171 }
 172
 173 void CreateHistos(TObjArray* harr){ // creates two histograms for each pattern in the DB
 174
 175   printf("\n\n Creating histograms...   ");
 176
 177   Int_t histonumber=0;
 178
 179   if(!harr) harr = &histoArr;
 180
 181   for(Int_t i=0; i<nPatterns; i++){
 182
 183     TH1* h = new TH1F(Form("hX%d", i),Form("#DeltaX for pattID %d", i),100,-30,30);
 184     h->SetDirectory(0);
 185     harr->AddAtAndExpand(h, i);
 186     histonumber++;
 187   }
 188
 189   for(Int_t i=0; i<nPatterns; i++){
 190
 191     TH1* h = new TH1F(Form("hZ%d", i),Form("#DeltaZ for pattID %d", i),100,-30,30);
 192     h->SetDirectory(0);
 193     harr->AddAtAndExpand(h, i+nPatterns);
 194     histonumber++;
 195   }
 196
 197   printf(" %d histograms created, corresponding to %d patterns\n\n", histonumber, histonumber/2);
 198 }
 199
 200 TH1* GetHistoX(Int_t pattID, TObjArray* harr){
 201
 202   TH1* h=0;
 203   if(!harr) harr = &histoArr;
 204
 205   h=(TH1*)harr->At(pattID);
 206   if (!h) {printf("Unknown histo id=%d\n",pattID); exit(1);}
 207   return h;
 208 }
 209
 210 TH1* GetHistoZ(Int_t pattID, TObjArray* harr){
 211
 212   Int_t zID=nPatterns+pattID;
 213
 214   TH1* h=0;
 215   if(!harr) harr = &histoArr;
 216
 217   h=(TH1*)harr->At(zID);
 218   if (!h) {printf("Unknown histo id=%d\n",zID); exit(1);}
 219   return h;
 220 }
 221
 222 void DrawReport(TObjArray* harr)
 223 {
 224   // plot all the nPatterns cluster
 225   for (int i=0;i<nPatterns;i++) {
 226     DrawNP(i,harr);
 227   }
 228 }
 229
 230 void DrawNP(int pattID, TObjArray* harr)
 231 {
 232   static TF1* gs = new TF1("gs","gaus",-500,500);
 233   if (!harr) harr = &histoArr;
 234
 235   printf("\n\nProcessing #DeltaX of pattern %d...",pattID);
 236
 237   Int_t fitStatusX=0;
 238
 239   TH1* hdx = (TH1*)harr->At(pattID);
 240   hdx->GetXaxis()->SetTitle("#DeltaX, #mum");
 241   hdx->SetTitle(Form("#DeltaX for pattID %d",pattID));
 242   gs->SetParameters(hdx->GetMaximum(),hdx->GetMean(),hdx->GetRMS());
 243   if((hdx->GetEntries())<100) fitStatusX = hdx->Fit("gs","Ql0");
 244   else fitStatusX = hdx->Fit("gs","Q0");
 245
 246   if(fitStatusX==0){
 247     Double_t px1 = gs->GetParameter(1);
 248     DeltaXmean[pattID]=px1;
 249     Double_t ex1 = gs->GetParError(1);
 250     DeltaXmeanErr[pattID]=ex1;
 251     Double_t px2 = gs->GetParameter(2);
 252     DeltaXsigma[pattID]=px2;
 253     Double_t ex2 = gs->GetParError(2);
 254     DeltaXsigmaErr[pattID]=ex2;
 255     Double_t ChiSqx = gs->GetChisquare();
 256     Chi2x[pattID] = ChiSqx;
 257     Int_t varNDFx = gs->GetNDF();
 258     if(varNDFx>=0)
 259     NDFx[pattID] = varNDFx;
 260     else{
 261       /*if(varNDFx==0){
 262         printf("\n\nNDF = 0, Chi2 = %e for pattern %d\n\n",ChiSqx,pattID);
 263         err3counter++;
 264       }*/
 265     NDFx[pattID]=-1;
 266     }
 267   }
 268   else{
 269     DeltaXmean[pattID]=0;
 270     DeltaXmeanErr[pattID]=0;
 271     DeltaXsigma[pattID]=0;
 272     DeltaXsigmaErr[pattID]=0;
 273     Chi2x[pattID] = -1;
 274   }
 275
 276   printf("done!!\n\n");
 277
 278   printf("Processing #DeltaZ of pattern %d... ",pattID);
 279
 280   Int_t fitStatusZ=0;
 281
 282   TH1* hdz = (TH1*)harr->At(pattID+nPatterns);
 283   hdz->SetTitle(Form("#DeltaZ for pattID %d",pattID));
 284   hdz->GetXaxis()->SetTitle("#DeltaZ, #mum");
 285   gs->SetParameters(hdz->GetMaximum(),hdz->GetMean(),hdz->GetRMS());
 286   if((hdz->GetEntries())<100) fitStatusZ = hdz->Fit("gs","Ql0");
 287   else fitStatusZ = hdz->Fit("gs","Q0");
 288
 289   if(fitStatusZ==0){
 290     Double_t pz1=gs->GetParameter(1);
 291     DeltaZmean[pattID]=pz1;
 292     Double_t ez1 = gs->GetParError(1);
 293     DeltaZmeanErr[pattID]=ez1;
 294     Double_t pz2 = gs->GetParameter(2);
 295     DeltaZsigma[pattID]=pz2;
 296     Double_t ez2 = gs->GetParError(2);
 297     DeltaZsigmaErr[pattID]=ez2;
 298     Double_t ChiSqz = gs->GetChisquare();
 299     Chi2z[pattID]=ChiSqz;
 300     Int_t varNDFz = gs->GetNDF();
 301     if(varNDFz>=0)
 302     NDFz[pattID] = varNDFz;
 303     else
 304     NDFz[pattID] = -1;
 305   }
 306   else{
 307     DeltaZmean[pattID]=0;
 308     DeltaZmeanErr[pattID]=0;
 309     DeltaZsigma[pattID]=0;
 310     DeltaZsigmaErr[pattID]=0;
 311     Chi2z[pattID] = -1;
 312   }
 313
 314   /*if((NDFz[pattID]!=NDFx[pattID])&&fitStatusZ!=fitStatusX){
 315     printf("\n\nERR: neither NDF nor FitStatus correspond!\n\n");
 316     err1counter++;
 317   }
 318   else if(NDFz[pattID]!=NDFx[pattID]){
 319     printf("\n\nERR: NDFx and NDFz donnot correspond!\n\n");
 320     err2counter++;
 321   }*/
 322
 323   printf("done!!\n\n");
 324
 325 }
 326
 327 void UpDateDB(const char* outDBFile){
 328
 329   printf("\n\n\nStoring data in the DataBase\n\n\n");
 330
 331   TString outDBFileName = outDBFile;
 332   if (outDBFileName.IsNull()) outDBFileName = "clusterTopology2.root";
 333
 334   TVectorF arrDeltaZmean(nPatterns);
 335   TVectorF arrDeltaZmeanErr(nPatterns);
 336   TVectorF arrDeltaXmean(nPatterns);
 337   TVectorF arrDeltaXmeanErr(nPatterns);
 338   TVectorF arrDeltaZsigma(nPatterns);
 339   TVectorF arrDeltaZsigmaErr(nPatterns);
 340   TVectorF arrDeltaXsigma(nPatterns);
 341   TVectorF arrDeltaXsigmaErr(nPatterns);
 342   TVectorF arrChi2x(nPatterns);
 343   TVectorF arrNDFx(nPatterns);
 344   TVectorF arrChi2z(nPatterns);
 345   TVectorF arrNDFz(nPatterns);
 346
 347   for(Int_t ID=0; ID<nPatterns; ID++){
 348
 349     printf("Processing pattern %d... ", ID);
 350
 351     arrDeltaZmean[ID]=DeltaZmean[ID];
 352     arrDeltaZmeanErr[ID]=DeltaZmeanErr[ID];
 353     arrDeltaXmean[ID]=DeltaXmean[ID];
 354     arrDeltaXmeanErr[ID]=DeltaXmeanErr[ID];
 355     arrDeltaZsigma[ID]=DeltaZsigma[ID];
 356     arrDeltaZsigmaErr[ID]=DeltaZsigmaErr[ID];
 357     arrDeltaXsigma[ID]=DeltaXsigma[ID];
 358     arrDeltaXsigmaErr[ID]=DeltaXsigmaErr[ID];
 359     arrChi2x[ID]=Chi2x[ID];
 360     arrNDFx[ID]=NDFx[ID];
 361     arrChi2z[ID]=Chi2z[ID];
 362     arrNDFz[ID]=NDFz[ID];
 363
 364     printf("done!!\n\n");
 365   }
 366
 367   TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
 368   flDB->WriteObject(pattDB,"TopDB","kSingleKey");
 369   flDB->WriteObject(NPix,"NPix","kSingleKey");
 370   flDB->WriteObject(NRow,"NRow","kSingleKey");
 371   flDB->WriteObject(NCol,"NCol","kSingleKey");
 372   flDB->WriteObject(pattFR, "TopFreq","kSingleKey");
 373   flDB->WriteObject(xCentrePix,"xCOG","kSingleKey");
 374   flDB->WriteObject(xCentreShift,"xShift","kSingleKey");
 375   flDB->WriteObject(zCentrePix,"zCOG","kSingleKey");
 376   flDB->WriteObject(zCentreShift,"zShift","kSingleKey");
 377   flDB->WriteObject(&arrDeltaZmean,"DeltaZmean","kSingleKey");
 378   flDB->WriteObject(&arrDeltaZmeanErr,"DeltaZmeanErr","kSingleKey");
 379   flDB->WriteObject(&arrDeltaXmean,"DeltaXmean","kSingleKey");
 380   flDB->WriteObject(&arrDeltaXmeanErr,"DeltaXmeanErr","kSingleKey");
 381   flDB->WriteObject(&arrDeltaZsigma,"DeltaZsigma","kSingleKey");
 382   flDB->WriteObject(&arrDeltaZsigmaErr,"DeltaZsigmaErr","kSingleKey");
 383   flDB->WriteObject(&arrDeltaXsigma,"DeltaXsigma","kSingleKey");
 384   flDB->WriteObject(&arrDeltaXsigmaErr,"DeltaXsigmaErr","kSingleKey");
 385   flDB->WriteObject(&arrChi2x,"Chi2x","kSingleKey");
 386   flDB->WriteObject(&arrChi2z,"Chi2z","kSingleKey");
 387   flDB->WriteObject(&arrNDFx,"NDFx","kSingleKey");
 388   flDB->WriteObject(&arrNDFz,"NDFz","kSingleKey");
 389   //
 390   flDB->Close();
 391   delete flDB;
 392
 393   printf("\n\nDB Complete!!\n\n");
 394 }