macros developed by Luca Barioglio for patterns analysis
[u/mrichter/AliRoot.git] / 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 }