macros developed by Luca Barioglio for patterns analysis
[u/mrichter/AliRoot.git] / ITS / UPGRADE / compression / PatternAnalysis / groupdef.C
1 /*
2 In this macro we create groups of patterns with similar characteristics. First we load the pattern database with information about MC-hit.
3 We do not group patterns with a frequency below a treshold (frequencyTreshold).
4 We define 7 different ID. 2 for the sigma of the MChit-COG, x and z direction (sigmaXID and sigma ZID);
5 2 for the distance between COG and centre of the pixel, x and z (shiftXID and shiftZID);
6 2 for the distance between MChit and centre of the pixel, x and z direction (biasXID ad biasZID).
7 1 ID for the number of fired pixels.
8 We define the binning (number/width of bins) and assign the previous IDs.
9 We decide which the grouping method we want to use: sigma, shift and number of pixels (kShift) or sigma, bias and number of pixels (kBias)
10 Finally we assign group ID. Frequent patterns (above the treshold), form a one-pattern group. The tohers are in the same group if they have
11 the same IDs.
12 A fil.txt is print qith he inormation for each apttern, including the IDs and the pattID of the patterns in the same group.
13
14 */
15
16 #if !defined(__CINT__) || defined(__MAKECINT__)
17 #include "TObjArray.h"
18 #include "TFile.h"
19 #include "TTree.h"
20 #include "TH1.h"
21 #include "TH2.h"
22 #include "TF1.h"
23 #include "TArrayI.h"
24 #include "TArrayF.h"
25 #include "TCanvas.h"
26 #include "TPad.h"
27 #include "TPavesText.h"
28 #include "TLatex.h"
29 #include "TBits.h"
30 #include "TGraph.h"
31 #include "TStopwatch.h"
32 #include "TMath.h"
33 #include "../ITS/UPGRADE/AliITSUClusterPix.h"
34 #include "../ITS/UPGRADE/AliITSURecoLayer.h"
35 #include "../ITS/UPGRADE/AliITSURecoDet.h"
36 #include "../ITS/UPGRADE/AliITSUHit.h"
37 #include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
38 #include "AliITSsegmentation.h"
39 #include "AliGeomManager.h"
40
41 #endif
42
43 TObjArray histoArr;
44 TObjArray *pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
45 TVectorF* pattFR=0; //frequency of the pattern in the database
46 TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
47 TVectorF* zCentrePix=0;
48 TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
49 TVectorF* zCentreShift=0;
50 TVectorF* NPix=0;//Number of fired pixels
51 TVectorF* NRow=0;//Number of rows of the pattern
52 TVectorF* NCol=0;//Number of columns of the pattern
53 TVectorF* DeltaZmean=0; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
54 TVectorF* DeltaXmean=0; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
55 TVectorF* DeltaZsigma=0; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
56 TVectorF* DeltaXsigma=0; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
57 TVectorF* DeltaZmeanErr=0;
58 TVectorF* DeltaXmeanErr=0;
59 TVectorF* DeltaZsigmaErr=0;
60 TVectorF* DeltaXsigmaErr=0;
61
62 //Defining the ID to create the groups.
63 TArrayI sigmaXID;
64 TArrayI sigmaZID;
65 TArrayI shiftXID;
66 TArrayI shiftZID;
67 TArrayI biasXID;
68 TArrayI biasZID;
69 TArrayI NPixID;
70 TArrayI groupID;
71 Float_t totalGroupFreq=0;
72
73 void LoadDB(const char* namefile);
74
75 Int_t nPatterns=0;
76
77 //Defining the frequency treshold under which ti group patterns
78 Float_t frequencyTreshold = 0.001;
79 //Defining the bins
80 Int_t NumberofSigmaXbins=1000;
81 Int_t NumberofSigmaZbins=1000;
82 Int_t NumberofShiftXbins=20;
83 Int_t NumberofShiftZbins=20;
84 Int_t NumberofBiasXbins=2000;
85 Int_t NumberofBiasZbins=2000;
86 //Defining the boundaries of the bins concerning the number of pixel
87 const Int_t     limitsnumber = 4;
88 //Defining the width of the bins to store patterns with similar sigma/shift/bias (in micron)
89 Float_t sigmaXwidth=3;
90 Float_t sigmaZwidth=3;
91 Float_t biasXwidth=5;
92 Float_t biasZwidth=5;
93 Float_t shiftXwidth=1./NumberofShiftXbins; //(fraction of the pitch)
94 Float_t shiftZwidth=1./NumberofShiftZbins; //(fraction of the pitch)
95
96 Int_t nPixlimits[limitsnumber] = {4,10,25,50};
97
98 enum{kShift=0, kBias=1};
99 //Select kShift to group wrt sigma & COG-centreOfThePixel distance
100 //Select kBias to group wrt sigma & MChit-centreOfThePixel distance
101
102 Int_t groupingMethod = kShift;
103
104 Int_t tempgroupID=0;
105
106 void groupdef(){
107
108         //Importing Data
109         LoadDB("clusterTopology2.root");
110
111         //Define moudule segmentation
112
113         AliGeomManager::LoadGeometry("geometry.root");
114         AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
115         AliITSUClusterPix::SetGeom(gm);
116         const AliITSsegmentation* segm = gm->GetSegmentation(0);
117         Float_t pitchX = segm->Dpx(0)*10000; // for pitch in X in micron
118         printf("pitchX: %f\n",pitchX);
119         Float_t pitchZ = segm->Dpz(0)*10000; // for pitch in Z in micron
120         printf("pitchX: %f\n",pitchX);
121         delete gm;
122
123         //Setting the number of patterns
124         nPatterns = (pattDB->GetEntriesFast());
125
126         sigmaXID.Set(nPatterns);
127         sigmaZID.Set(nPatterns);
128
129         //Assign -2 to frequent patterns, not to group, and -1 to rare clusters, -3 to patterns with sigma set to zero
130         for(Int_t ID=0; ID<nPatterns; ID++){
131                 printf("Assign temporary sigma ID to pattern %d... ",ID);
132                 if((*pattFR)[ID]>frequencyTreshold){
133                         sigmaXID[ID]=-2;//In order not to consider it within the groups
134                         sigmaZID[ID]=-2;//In order not to consider it within the groups
135                 }
136                 else{
137                         sigmaXID[ID]=-1;
138                         sigmaZID[ID]=-1;
139                         totalGroupFreq+=(*pattFR)[ID];
140                 }
141                 printf("done\n\n");
142         }
143
144         //Assign to similar patterns the same sigmaXID
145         for(Int_t i=0; i<nPatterns; i++){
146                 printf("Assign sigmaXID to pattern %d... ",i);
147                 if(sigmaXID[i]==-1){
148                         if((*DeltaXsigma)[i]==0){
149                                 sigmaXID[i]=-3; // In order not to cosider it within the groups
150                         }
151                         else{
152                                 for(Int_t j=0; j<NumberofSigmaXbins; j++){
153                                         if(j*sigmaXwidth < (*DeltaXsigma)[i] && (*DeltaXsigma)[i]<= (j+1)*sigmaXwidth){
154                                                 sigmaXID[i]=j+1;
155                                                 break;
156                                         }
157                                 }
158                         }
159                 }
160                 printf("done!!\n\n");
161         }
162
163         //Assign to similar patterns the same sigmaZID
164         for(Int_t i=0; i<nPatterns; i++){
165                 printf("Assign sigmaZID to pattern %d... ",i);
166                 if(sigmaZID[i]==-1){
167                         if((*DeltaZsigma)[i]==0){
168                                 sigmaZID[i]=-3; // In order not to cosider it within the groups
169                         }
170                         else{
171                                 for(int j=0; j<NumberofSigmaZbins ; j++){
172                                         if(j*sigmaZwidth < (*DeltaZsigma)[i] && (*DeltaZsigma)[i]<= (j+1)*sigmaZwidth){
173                                                 sigmaZID[i]=j+1;
174                                                 break;
175                                         }
176                                 }
177                         }
178                 }
179                 printf("done!!\n\n");
180         }
181
182         //assigning shiftID
183
184         shiftXID.Set(nPatterns);
185         shiftZID.Set(nPatterns);
186
187         for(Int_t i=0; i<nPatterns; i++){
188                 printf("Assign shiftXID to pattern %d... ",i);
189                 
190                 for(int j=-NumberofShiftXbins/2; j<NumberofShiftXbins/2; j++){
191                         
192                         if(j*shiftXwidth < (*xCentreShift)[i] && (*xCentreShift)[i]<= (j+1)*shiftXwidth){
193                                 shiftXID[i]=j+1;
194                                 printf("done!!\n\n");
195                                 break;
196                         }       
197                 }
198         }
199
200         for(Int_t i=0; i<nPatterns; i++){
201                 printf("Assign shiftZID to pattern %d... ",i);
202                 
203                 for(int j=-NumberofShiftZbins/2; j<NumberofShiftZbins/2; j++){
204                         if(j*shiftZwidth < (*zCentreShift)[i] && (*zCentreShift)[i]<= (j+1)*shiftZwidth){
205                                 shiftZID[i]=j+1;
206                                 printf("done!!\n\n");
207                                 break;
208                         }
209                 }       
210         }
211         
212         //assigning BiasID
213
214         biasXID.Set(nPatterns);
215         biasZID.Set(nPatterns);
216
217         //Setting all the bias ID to zero
218
219         for(Int_t i=0; i<nPatterns; i++){
220                 biasXID[i]=0;
221                 biasZID[i]=0;
222         }
223
224         for(Int_t i=0; i<nPatterns; i++){
225                 printf("Assign biasXID to pattern %d... ",i);
226                 for(Int_t j=-NumberofBiasXbins/2; j<NumberofBiasXbins/2; j++){
227                         if(j*biasXwidth < ((*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX))
228                                 && ((*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX))<= (j+1)*biasXwidth){
229                         biasXID[i]=j+1;
230                         break;
231                         }
232                 }
233                 printf("done!!\n\n");
234         }
235
236         for(Int_t i=0; i<nPatterns; i++){
237                 biasXID[i]=0;
238                 biasZID[i]=0;
239         }
240
241         for(Int_t i=0; i<nPatterns; i++){
242                 printf("Assign biasZID to pattern %d... ",i);
243                 for(Int_t j=-NumberofBiasZbins/2; j<NumberofBiasZbins/2; j++){
244                         if(j*biasZwidth < ((*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ)) 
245                                 && ((*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ))<= (j+1)*biasZwidth){
246                         biasZID[i]=j+1;
247                         break;
248                         }
249                 }
250                 printf("done!!\n\n");
251         }
252
253         
254
255         //Assigning NPixID
256
257         NPixID.Set(nPatterns);
258
259         for(Int_t i=0; i<nPatterns; i++){
260                 printf("Assigning NPixID to pattern %d...", i);
261                 if((*NPix)[i]<=nPixlimits[0]){
262                         NPixID[i]=0;
263                         printf("done!!\n\n");
264                 }
265                 else if(nPixlimits[0]<(*NPix)[i] && (*NPix)[i]<=nPixlimits[limitsnumber-1]){
266                         for(Int_t j=0; j<nPatterns; j++ ){
267                                 if(nPixlimits[j]<(*NPix)[i] && (*NPix)[i]<=nPixlimits[j+1]){
268                                         NPixID[i]=j+1;
269                                         printf("done!!\n\n");
270                                         break;
271                                 }
272                         }
273                 }
274                 else if((*NPix)[i]>nPixlimits[limitsnumber-1]){
275                         NPixID[i]=limitsnumber+1;
276                         printf("done!!\n\n");
277                 }
278         }
279
280         //Assigning groupID
281
282         groupID.Set(nPatterns);
283
284         //Assign -2 to frequent patterns, not to group, and -1 to rare clusters
285         for(Int_t ID=0; ID<nPatterns; ID++){
286                 printf("Assign temporary group ID to pattern %d... ",ID);
287                 groupID[ID]=-1;
288                 printf("done\n\n");
289         }
290         Int_t k=0;
291         while((*pattFR)[k]>frequencyTreshold){
292                 groupID[k]=tempgroupID;
293                 tempgroupID++;
294                 k++;
295         }
296
297         if(groupingMethod==kShift){
298                 for(Int_t i=0; i<nPatterns; i++){
299                         if(groupID[i]!=-1) continue;    
300                         groupID[i]=tempgroupID;
301                         printf("Assigning group ID %d... ",tempgroupID);
302                         for(Int_t j=i+1; j<nPatterns; j++){
303                                 if(sigmaXID[j]==-3){
304                                         groupID[j]=-1;
305                                         continue;
306                                 }
307                                 else if(sigmaXID[j]==sigmaXID[i] && sigmaZID[j]==sigmaZID[i] && 
308                                         shiftXID[j]==shiftXID[i] && shiftZID[j]==shiftZID[i] &&
309                                         NPixID[i]==NPixID[j]) groupID[j]=tempgroupID;
310                         }
311                         printf("done!!\n\n");
312                         tempgroupID++;
313                 }
314         }
315         else if(groupingMethod==kBias){
316                 for(Int_t i=0; i<nPatterns; i++){
317                         if(groupID[i]!=-1) continue;
318                         groupID[i]=tempgroupID;
319                         printf("Assigning group ID %d... ",tempgroupID);
320                         for(Int_t j=i+1; j<nPatterns; j++){
321                                 if(sigmaZID[j]==-3){
322                                         groupID[j]=-1;
323                                         continue;
324                                 }
325                                 else if(sigmaXID[j]==sigmaXID[i] && sigmaZID[j]==sigmaZID[i]
326                                         && biasXID[j]==biasXID[i] && biasZID[j]==biasZID[i]
327                                         && NPixID[i]==NPixID[j]) groupID[j]=tempgroupID;
328                         }
329                         printf("done!!\n\n");
330                         tempgroupID++;
331                 }
332         }
333
334         ofstream a("groupdef.txt");
335                 
336         //setw(55) << "patterns in the group\n" << endl;
337
338         a << Form("A NEGATIVE ID means that the pattern is not in a group.") << endl << endl <<
339                 Form("EXCEPTION: biasID CAN be negative") << endl <<
340                 "\n\n......................................................................................." << 
341                 "................................................................................................\n\n";
342
343         for(int i=0; i<nPatterns; i++){
344
345                 printf("Writing info about pattern %d ...", i);
346
347                 a <<  setw(30) << Form("pattID: %d",i) <<  setw(30) << Form("freq: %f", (*pattFR)[i]) << setw(30) << Form("NPix: %d",Int_t((*NPix)[i]))<< setw(45) << Form("NRow: %d ",Int_t((*NRow)[i])) << setw(45) <<
348                 Form("NCol: %d", Int_t((*NCol)[i])) << endl << endl
349                 << setw(45) << Form("DeltaXmean: %f (%f)",(*DeltaXmean)[i],(*DeltaXmeanErr)[i]) << setw(45) << 
350                 Form("DeltaZmean: %f (%f)",(*DeltaZmean)[i],(*DeltaZmeanErr)[i])<<
351                 setw(45) << Form("DeltaXsigma: %f (%f)",(*DeltaXsigma)[i],(*DeltaXsigmaErr)[i]) << setw(45) <<
352                 Form("DeltaZsigma: %f (%f)",(*DeltaZsigma)[i],(*DeltaZsigmaErr)[i]) << endl << endl << setw(30) << 
353                 Form("xShift: %f",(*xCentreShift)[i]) <<  setw(45) << Form("zShift: %f",(*zCentreShift)[i]) << setw(45) << Form("BiasX(MChit-centrePix): %f",(*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX)) << setw(45) <<
354                 Form("BiasZ(MChit-centrePix): %f",(*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ)) << endl << endl <<
355                 setw(30) << Form("sigmaXID: %d",sigmaXID[i]) << setw(15) << Form("sigmaZID: %d",sigmaZID[i]) << setw(15) <<
356                 Form("shiftXID: %d",shiftXID[i]) << setw(15) << Form("shiftZID: %d",shiftZID[i]) << setw(15)<< Form("biasXID: %d",biasXID[i]) <<
357                 setw(15) << Form("biasZID: %d",biasZID[i]) << setw(15) << Form("NPixID: %d", NPixID[i]) << endl << endl << setw(30) << Form("groupID: %d", groupID[i])
358                 << endl << endl << setw(30) << "patterns in this group:\n\n";
359
360                 Int_t matchcounter=0;
361
362                 for(int j=0; j<nPatterns;  j++){
363                         if(matchcounter!=0 && (matchcounter%20)==0){
364                                 a << endl;
365                                 matchcounter=0;
366                         }
367                         if(groupID[j]==groupID[i]) {
368                                 a << j << ", ";
369                                 matchcounter++;
370                         }
371                 }
372                 
373                 a << "\n\n......................................................................................." << 
374                 "................................................................................................\n\n";
375
376
377                 printf("done!\n\n");
378         }
379
380         a.close();
381
382         printf("%d groups found!!!!\n\n",tempgroupID);
383
384         printf("\n\nThe total frequency of the patterns in group is %f\n\n",totalGroupFreq);
385
386         /*
387
388         TVectorF groupDeltaX(tempGID+100);
389         TVectorF groupShiftX(tempGID+100);
390         TVectorF groupDeltaZ(tempGID+100);
391         TVectorF groupShiftZ(tempGID+100);
392         TVectorF patternNum(tempGID+100);
393         groupDeltaX.Zero();
394         groupShiftX.Zero();
395         groupDeltaZ.Zero();
396         groupShiftZ.Zero();
397         patternNum.Zero();
398
399         ofstream b("groupDB.txt");
400
401         b<<setw(15)<<"groupID"<<setw(25)<<"patterns in the group"<<setw(15)<<
402         "DeltaX"<<setw(15)<<"ShiftX"<<setw(15)<<"DeltaZ"<<setw(15)<<"ShiftZ\n"<<endl;
403
404         TCanvas* c = new TCanvas("c","Patterns groups",900,600);
405         TH1F* h = new TH1F("h","Patterns groups",tempGID,-0.5, tempGID-0.5);
406         h->GetXaxis()->SetTitle("group ID");
407         h->SetStats(0);
408
409         for(Int_t gid=0; gid<tempGID; gid++){
410
411                 Int_t PatternNumberInGroup=0;
412                 Float_t freqSum=0.;//It is the sum of the frequencies of the patterns in the same group
413
414                 for(Int_t pattID=0; pattID<nPatterns; pattID++){
415                         if (groupID[pattID]==gid)
416                         {
417                                 groupDeltaX[gid]+=(*DeltaXsigma)[pattID]*(*pattFR)[pattID];
418                                 groupShiftX[gid]+=(*xCentreShift)[pattID]*(*pattFR)[pattID];
419                                 groupDeltaZ[gid]+=(*DeltaZsigma)[pattID]*(*pattFR)[pattID];
420                                 groupShiftZ[gid]+=(*zCentreShift)[pattID]*(*pattFR)[pattID];
421                                 freqSum+=(*pattFR)[pattID];
422                                 PatternNumberInGroup++;
423
424                                 h->Fill(gid);
425                         }
426                 }
427
428                 groupDeltaX[gid]=groupDeltaX[gid]/freqSum;
429                 groupShiftX[gid]=groupShiftX[gid]/freqSum;
430                 groupDeltaZ[gid]=groupDeltaZ[gid]/freqSum;
431                 groupShiftZ[gid]=groupShiftZ[gid]/freqSum;
432                 patternNum[gid]=PatternNumberInGroup;
433
434                 b<<setw(15)<<gid<<setw(25)<<patternNum[gid]<<setw(15)<<
435                 groupDeltaX[gid]<<setw(15)<<groupShiftX[gid]<<setw(15)<<
436                 groupDeltaZ[gid]<<setw(15)<<groupShiftZ[gid]<<"\n"<<endl;
437         }
438
439         c->cd();
440         h->Draw();
441
442         TPavesText* info = new TPavesText(0.5,1,0.5,1,1,"nb");
443         info->AddText(Form("Number of groups: %d",tempGID));
444         info->AddText("#delta_{X}<5 & #delta_{Z}<5");
445         info->AddText("#DeltaX<0.05 & #DeltaZ<0.05");
446         c->cd();
447         info->Draw();
448
449         b.close();
450         */
451 }
452
453 void LoadDB(const char* fname){
454
455   printf("\n\nLoading DB... ");
456   // load database
457   TFile* fl = TFile::Open(fname);
458   if(!fl){printf("Could not find %s",fname); exit(1);}
459   pattDB = (TObjArray*)fl->Get("TopDB");
460   pattFR = (TVectorF*)fl->Get("TopFreq");
461   xCentrePix =(TVectorF*)fl->Get("xCOG");
462   zCentrePix =(TVectorF*)fl->Get("zCOG");
463   xCentreShift =(TVectorF*)fl->Get("xShift");
464   zCentreShift =(TVectorF*)fl->Get("zShift");
465   NPix =(TVectorF*)fl->Get("NPix");
466   NCol =(TVectorF*)fl->Get("NCol");
467   NRow =(TVectorF*)fl->Get("NRow");
468   DeltaXmean =(TVectorF*)fl->Get("DeltaXmean");
469   DeltaZmean =(TVectorF*)fl->Get("DeltaZmean");
470   DeltaXsigma =(TVectorF*)fl->Get("DeltaXsigma");
471   DeltaZsigma =(TVectorF*)fl->Get("DeltaZsigma");
472   DeltaXmeanErr =(TVectorF*)fl->Get("DeltaXmeanErr");
473   DeltaZmeanErr =(TVectorF*)fl->Get("DeltaZmeanErr");
474   DeltaZsigmaErr =(TVectorF*)fl->Get("DeltaXsigmaErr");
475   DeltaXsigmaErr =(TVectorF*)fl->Get("DeltaZsigmaErr");
476   printf("done!!\n\n");
477 }