Reduced QA output (Yves)
[u/mrichter/AliRoot.git] / ITS / AliITSOnlineSPDphysAnalyzer.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ////////////////////////////////////////////////////////////
17 // Author: Henrik Tydesjo                                 //
18 // This class is used in the detector algorithm framework //
19 // to process the data stored in special container files  //
20 // (see AliITSOnlineSPDphys).                             //
21 ////////////////////////////////////////////////////////////
22
23 #include "AliITSOnlineSPDphysAnalyzer.h"
24 #include "AliITSOnlineSPDphys.h"
25 #include "AliITSOnlineCalibrationSPDhandler.h"
26 #include "AliITSRawStreamSPD.h"
27 #include "AliITSIntMap.h"
28 #include <TStyle.h>
29 #include <TMath.h>
30 #include <TF1.h>
31 #include <TGraph.h>
32 #include <TH2F.h>
33 #include <TError.h>
34 #include <iostream>
35 #include <fstream>
36
37 AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(const Char_t *fileName, AliITSOnlineCalibrationSPDhandler* handler, Bool_t readFromGridFile) :
38   fFileName(fileName),fPhysObj(NULL),fHandler(handler),
39   fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
40   fNrEqHits(0),fbDeadProcessed(kFALSE),
41   fThreshNoisy(1e-9),fThreshDead(1e-9),
42   fMinEventsForNoisy(10000),fMinEventsForDead(10000),
43   fDefinitelyNoisyRatio(0.3),
44   fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)
45 {
46   // constructor  
47   Init(readFromGridFile);
48 }
49
50 AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(AliITSOnlineSPDphys* physObj, AliITSOnlineCalibrationSPDhandler* handler) :
51   fFileName("test.root"),fPhysObj(NULL),fHandler(handler),
52   fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
53   fNrEqHits(0),fbDeadProcessed(kFALSE),
54   fThreshNoisy(1e-9),fThreshDead(1e-9),
55   fMinEventsForNoisy(10000),fMinEventsForDead(10000),
56   fDefinitelyNoisyRatio(0.3),
57   fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)
58 {
59   // alt constructor
60   fPhysObj = physObj;
61 }
62
63 AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(const AliITSOnlineSPDphysAnalyzer& handle) :
64   fFileName("test.root"),fPhysObj(NULL),fHandler(NULL),
65   fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
66   fNrEqHits(0),fbDeadProcessed(kFALSE),
67   fThreshNoisy(1e-9),fThreshDead(1e-9),
68   fMinEventsForNoisy(10000),fMinEventsForDead(10000),
69   fDefinitelyNoisyRatio(0.3),
70   fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)
71 {
72   // copy constructor, only copies the filename and params (not the processed data)
73   fFileName=handle.fFileName;
74   fThreshNoisy = handle.fThreshNoisy;
75   fThreshDead = handle.fThreshDead;
76   fMinEventsForNoisy = handle.fMinEventsForNoisy;
77   fMinEventsForDead = handle.fMinEventsForDead;
78   fDefinitelyNoisyRatio = handle.fDefinitelyNoisyRatio;
79   fMinNrEqHitsForDeadChips = handle.fMinNrEqHitsForDeadChips;
80   fRatioToMeanForInefficientChip = handle.fRatioToMeanForInefficientChip;
81
82   Init();
83 }
84
85 AliITSOnlineSPDphysAnalyzer::~AliITSOnlineSPDphysAnalyzer() {
86   // destructor
87   if (fPhysObj!=NULL) delete fPhysObj;
88 }
89
90 AliITSOnlineSPDphysAnalyzer& AliITSOnlineSPDphysAnalyzer::operator=(const AliITSOnlineSPDphysAnalyzer& handle) {
91   // assignment operator, only copies the filename and params (not the processed data)
92   if (this!=&handle) {
93     if (fPhysObj!=NULL) delete fPhysObj;
94     
95     fFileName=handle.fFileName;
96     fThreshNoisy = handle.fThreshNoisy;
97     fThreshDead = handle.fThreshDead;
98     fMinEventsForNoisy = handle.fMinEventsForNoisy;
99     fMinEventsForDead = handle.fMinEventsForDead;
100     fDefinitelyNoisyRatio = handle.fDefinitelyNoisyRatio;
101     fMinNrEqHitsForDeadChips = handle.fMinNrEqHitsForDeadChips;
102     fRatioToMeanForInefficientChip = handle.fRatioToMeanForInefficientChip;
103
104     fPhysObj = NULL;
105     fHandler = NULL;
106     fNrEnoughStatChips = 0;
107     fNrDeadChips = 0;
108     fNrInefficientChips = 0;
109     fNrEqHits = 0;
110     fbDeadProcessed = kFALSE;
111
112     Init();    
113   }
114   return *this;
115 }
116
117 void AliITSOnlineSPDphysAnalyzer::Init(Bool_t readFromGridFile) {
118   // initialize container obj
119   if (!readFromGridFile) {
120     FILE* fp0 = fopen(fFileName.Data(), "r");
121     if (fp0 == NULL) {
122       return;
123     }
124     else {
125       fclose(fp0);
126     }
127   }
128   fPhysObj = new AliITSOnlineSPDphys(fFileName.Data(), readFromGridFile);
129 }
130
131 void AliITSOnlineSPDphysAnalyzer::SetParam(const Char_t *pname, const Char_t *pval) {
132   // set a parameter
133   TString name = pname;
134   TString val = pval;
135   //  printf("Setting Param %s  to %s\n",name.Data(),val.Data());
136   if (name.CompareTo("MistakeProbabilityNoisy")==0) {
137     Double_t mistakeProbabilityNoisy = val.Atof();
138     fThreshNoisy = mistakeProbabilityNoisy/(20*6*10*32*256);
139   }
140   else if (name.CompareTo("MistakeProbabilityDead")==0) {
141     Double_t mistakeProbabilityDead = val.Atof();
142     fThreshDead = mistakeProbabilityDead/(20*6*10*32*256);
143   }
144   else if (name.CompareTo("fMinEventsForNoisy")==0) {
145     fMinEventsForNoisy = val.Atoi();
146   }
147   else if (name.CompareTo("fMinEventsForDead")==0) {
148     fMinEventsForDead = val.Atoi();
149   }
150   else if (name.CompareTo("fDefinitelyNoisyRatio")==0) {
151     fDefinitelyNoisyRatio = val.Atof();
152   }
153   else if (name.CompareTo("fMinNrEqHitsForDeadChips")==0) {
154     fMinNrEqHitsForDeadChips = val.Atof();
155   }
156   else if (name.CompareTo("fRatioToMeanForInefficientChip")==0) {
157     fRatioToMeanForInefficientChip = val.Atof();
158   }
159   else {
160     Error("AliITSOnlineSPDphysAnalyzer::SetParam","Parameter %s in configuration file unknown.",name.Data());
161   }
162 }
163
164 void AliITSOnlineSPDphysAnalyzer::ReadParamsFromLocation(const Char_t *dirName) {
165   // opens file (default name) in dir dirName and reads parameters from it
166   TString paramsFileName = Form("%s/physics_params.txt",dirName);
167   ifstream paramsFile;
168   paramsFile.open(paramsFileName, ifstream::in);
169   if (paramsFile.fail()) {
170     printf("No config file (%s) present. Using default tuning parameters.\n",paramsFileName.Data());
171   }
172   else {
173     while(1) {
174       Char_t paramN[50];
175       Char_t paramV[50];
176       paramsFile >> paramN;
177       if (paramsFile.eof()) break;
178       paramsFile >> paramV;
179       SetParam(paramN,paramV);
180       if (paramsFile.eof()) break;
181     }
182     paramsFile.close();
183   }
184 }
185
186
187 UInt_t AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels() {
188   // process noisy pixel data , returns number of chips with enough statistics
189   if (fPhysObj==NULL) {
190     Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","No data!");
191     return 0;
192   }
193   // do we have enough events to even try the algorithm?
194   if (GetNrEvents() < fMinEventsForNoisy) {
195     Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Nr events (%d) < fMinEventsForNoisy (%d)!",GetNrEvents(),fMinEventsForNoisy);
196     return 0;
197   }
198   // handler should be initialized
199   if (fHandler==NULL) {
200     Error("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Calibration handler is not initialized!");
201     return 0;
202   }
203   
204   UInt_t nrEnoughStatChips = 0;
205
206   for (UInt_t hs=0; hs<6; hs++) {
207     for (UInt_t chip=0; chip<10; chip++) {
208
209       UInt_t nrPixels = 0;
210       UInt_t nrChipHits = 0;
211       UInt_t nrMostHits = 0;
212       for (UInt_t col=0; col<32; col++) {
213         for (UInt_t row=0; row<256; row++) {
214           UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
215           nrChipHits += nrHits;
216           //      if (nrHits>0) nrPixels++; // don't include pixels that might be dead
217           nrPixels++;
218           if (nrHits>fDefinitelyNoisyRatio*GetNrEvents()) {
219             fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
220             nrPixels--;
221             nrChipHits-=nrHits;
222           }
223           else {
224             if (nrMostHits<nrHits) nrMostHits=nrHits;
225           }
226         }
227       }
228
229       if (nrChipHits>0) { // otherwise there are for sure no noisy
230         // Binomial with n events and probability p for pixel hit
231         UInt_t n = GetNrEvents();
232         if (nrPixels>0 && n>0) {
233
234           Double_t p = (Double_t)nrChipHits/nrPixels/n;
235
236           // Bin(n,k=0):
237           Double_t bin = pow((Double_t)(1-p),(Double_t)n);
238           // Bin(n,k)
239           UInt_t k=1;
240           while ((bin>fThreshNoisy || k<n*p) && k<=n) {
241             k++;
242             bin = bin*(n-k+1)/k*p/(1-p);
243           }
244           
245           // can we find noisy pixels...?
246           if (k<=n) {
247             //      printf("eq %d , hs %d , chip %d : Noisy level = %d\n",GetEqNr(),hs,chip,k);
248             nrEnoughStatChips++;
249             // add noisy pixels to handler
250             UInt_t noiseLimit=k;
251             if (nrMostHits>=noiseLimit) {
252               for (UInt_t col=0; col<32; col++) {
253                 for (UInt_t row=0; row<256; row++) {
254                   UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
255                   if (nrHits >= noiseLimit) {
256                     fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
257                   }
258                 }
259               }
260             }
261           }
262         }
263
264       }
265
266     } // for chip
267   } // for hs
268
269   return nrEnoughStatChips;
270 }
271
272
273 UInt_t AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels() {
274   // process dead pixel data , returns number of chips with enough statistics
275   if (fPhysObj==NULL) {
276     Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","No data!");
277     return 0;
278   }
279
280   // do we have enough events to even try the algorithm?
281   if (GetNrEvents() < fMinEventsForDead) {
282     Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Nr events (%d) < fMinEventsForDead (%d)!",GetNrEvents(),fMinEventsForDead);
283     return 0;
284   }
285   // handler should be initialized
286   if (fHandler==NULL) {
287     Error("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Calibration handler is not initialized!");
288     return 0;
289   }
290
291   AliITSIntMap* possiblyDead  = new AliITSIntMap();
292   AliITSIntMap* possiblyIneff = new AliITSIntMap();
293
294   fNrEnoughStatChips = 0;
295   fNrDeadChips = 0;
296   fNrInefficientChips = 0;
297   UInt_t nrPossiblyDeadChips = 0;
298   fNrEqHits = 0;
299
300
301   for (UInt_t hs=0; hs<6; hs++) {
302     if (!fHandler->IsActiveHS(GetEqNr(),hs)) {
303       fNrDeadChips+=10;
304     }
305     else {
306       for (UInt_t chip=0; chip<10; chip++) {
307         if (!fHandler->IsActiveChip(GetEqNr(),hs,chip)) {
308           fNrDeadChips++;
309         }
310         else {
311           // perform search for individual dead pixels...
312           Bool_t good=kFALSE;
313
314           UInt_t nrPossiblyDeadPixels = 0;
315           UInt_t nrPixels = 0;
316           UInt_t nrChipHits = 0;
317           for (UInt_t col=0; col<32; col++) {
318             for (UInt_t row=0; row<256; row++) {
319               UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
320               nrChipHits += nrHits;
321               if (!fHandler->IsPixelNoisy(GetEqNr(),hs,chip,col,row)) {
322                 // don't include noisy pixels
323                 nrPixels++;
324                 if (nrHits==0) {
325                   nrPossiblyDeadPixels++;
326                 }
327                 else {
328                   fHandler->UnSetDeadPixel(GetEqNr(),hs,chip,col,row); // unset (no action unless dead before)
329                 }
330               }
331               else {
332                 nrChipHits -= nrHits; // this is needed when running offline (online nrHits should be 0 already)
333               }
334             }
335           }
336           fNrEqHits+=nrChipHits;
337
338           if (nrChipHits>0) {
339             // make sure the chip is not flagged as dead
340             fHandler->SetDeadChip(GetEqNr(),hs,chip,kFALSE);
341           }
342
343           if (nrPossiblyDeadPixels==0) {
344             // no need to see if we have enough statistics...
345             fNrEnoughStatChips++;
346             good=kTRUE;
347             //  printf("%3d",good);
348             //  if (chip==9) printf("\n");
349             continue;
350           }
351
352           if (nrChipHits==0) {
353             nrPossiblyDeadChips++;
354             possiblyDead->Insert(hs,chip);
355             good=kFALSE;
356             //  printf("%3d",good);
357             //  if (chip==9) printf("\n");
358             continue;
359           }
360
361           // Binomial with n events and probability p for pixel hit
362           UInt_t n = GetNrEvents();
363           if (nrPixels>0 && n>0) {
364
365             Double_t p = (Double_t)nrChipHits/nrPixels/n;
366
367             // probability of falsely assigning a dead pixel
368             Double_t falselyDeadProb = pow((Double_t)(1-p),(Double_t)n);
369             //      printf("falselyprob=%e\n",falselyDeadProb);
370
371             // can we find dead pixels...?
372             if (falselyDeadProb<fThreshDead) {
373               fNrEnoughStatChips++;
374               good=kTRUE;
375               // add dead pixels to handler
376               for (UInt_t col=0; col<32; col++) {
377                 for (UInt_t row=0; row<256; row++) {
378                   UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
379                   if (nrHits==0) {
380                     if (!fHandler->IsPixelNoisy(GetEqNr(),hs,chip,col,row)) {
381                       // don't include noisy pixels
382                       fHandler->SetDeadPixel(GetEqNr(),hs,chip,col,row);
383                     }
384                   }
385                 }
386               }
387             }
388             if (!good) {
389               // this might be an inefficient chip
390               possiblyIneff->Insert(hs*10+chip,nrChipHits);
391             }
392
393           }
394           else {
395             if (n>0) {
396               // this is a completely noisy chip... put in category enough stat
397               fNrEnoughStatChips++;
398               good=kTRUE;
399             }
400           }
401
402           //      printf("%3d",good);
403           //      if (chip==9) printf("\n");
404
405         }
406       } // for chip
407     }
408   } // for hs
409  
410
411   Int_t key,val;
412
413   // dead chips?
414   if (fNrEqHits>fMinNrEqHitsForDeadChips) {
415     while (possiblyDead->Pop(key,val)) {
416       fHandler->SetDeadChip(GetEqNr(),key,val,kFALSE);
417     }
418     fNrDeadChips+=nrPossiblyDeadChips;
419   }
420   delete possiblyDead;
421
422   // inefficient chips?
423   while (possiblyIneff->Pop(key,val)) {
424     if (val<fNrEqHits/60*fRatioToMeanForInefficientChip) {
425       fNrInefficientChips++;
426     }
427   }
428   delete possiblyIneff;
429
430
431   fbDeadProcessed = kTRUE;
432
433   return fNrEnoughStatChips;
434 }
435
436
437 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEnoughStatChips() {
438   // returns nr of enough stat chips
439   if (!fbDeadProcessed) ProcessDeadPixels();
440   return fNrEnoughStatChips;
441 }
442 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrDeadChips() {
443   // returns nr of dead chips
444   if (!fbDeadProcessed) ProcessDeadPixels();
445   return fNrDeadChips;
446 }
447 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrInefficientChips() {
448   // returns nr of inefficient chips
449   if (!fbDeadProcessed) ProcessDeadPixels();
450   return fNrInefficientChips;
451 }
452 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrNeedsMoreStatChips() {
453   // returns nr of needs more stat chips
454   if (!fbDeadProcessed) ProcessDeadPixels();
455   return 60-fNrEnoughStatChips-fNrDeadChips-fNrInefficientChips;
456 }
457
458 UInt_t AliITSOnlineSPDphysAnalyzer::GetEqNr() const {
459   // returns the eq nr of phys obj
460   if (fPhysObj!=NULL) return fPhysObj->GetEqNr(); 
461   else return 999;
462 }
463
464 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEvents() const {
465   // returns the nr of events of phys obj
466   if (fPhysObj!=NULL) return fPhysObj->GetNrEvents(); 
467   else return 0;
468 }
469
470 void AliITSOnlineSPDphysAnalyzer::Exponent(Double_t &val, Int_t &valExp) const {
471   // put double in format with val and exp so that 1<val<10 - The actual value is val*10e(valExp)
472   while (val>10) {
473     val/=10;
474     valExp++;
475   }
476   while (val<1) {
477     val*=10;
478     valExp--;
479   }
480 }
481
482
483
484 TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapTot() {
485   // creates and returns a pointer to a hitmap histo (half sector style a la spdmood)
486   if (fPhysObj==NULL) {
487     Error("AliITSOnlineSPDphysAnalyzer::GetHitMapTot","No data!");
488     return NULL;
489   }
490   TString histoname = Form("Eq %d",GetEqNr());
491   TH2F* fHitMapTot = new TH2F(histoname.Data(),histoname.Data(),32*10,-0.5,32*10-0.5,256*6,-0.5,256*6-0.5);
492   fHitMapTot->SetNdivisions(-10,"X");
493   fHitMapTot->SetNdivisions(-006,"Y");
494   fHitMapTot->SetTickLength(0,"X");
495   fHitMapTot->SetTickLength(0,"Y");
496   fHitMapTot->GetXaxis()->SetLabelColor(gStyle->GetCanvasColor());
497   fHitMapTot->GetYaxis()->SetLabelColor(gStyle->GetCanvasColor());
498   for (UInt_t hs=0; hs<6; hs++) {
499     for (UInt_t chipNr=0; chipNr<10; chipNr++) {
500       for (UInt_t col=0; col<32; col++) {
501         for (UInt_t row=0; row<256; row++) {
502           fHitMapTot->Fill(chipNr*32+col,(5-hs)*256+row,fPhysObj->GetHits(hs,chipNr,col,row));
503         }
504       }
505     }
506   }
507   return fHitMapTot;
508 }
509
510 TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapChip(UInt_t hs, UInt_t chip) {
511   // creates and returns a pointer to a hitmap histo (chip style a la spdmood)
512   if (fPhysObj==NULL) {
513     Error("AliITSOnlineSPDphysAnalyzer::GetHitMapChip","No data!");
514     return NULL;
515   }
516
517   TString histoName;
518   TString histoTitle;
519   histoName = Form("fChipHisto_%d_%d_%d", GetEqNr(), hs, chip);
520   histoTitle = Form("Eq ID %d, Half Stave %d, Chip %d", GetEqNr(), hs, chip);
521
522   TH2F *returnHisto = new TH2F(histoName.Data(), histoTitle.Data(), 32, -0.5, 31.5, 256, -0.5, 255.5);
523   returnHisto->SetMinimum(0);
524   for (UInt_t col=0; col<32; col++) {
525     for (UInt_t row=0; row<256; row++) {
526       returnHisto->Fill(col,row,fPhysObj->GetHits(hs,chip,col,row));
527     }
528   }
529
530   return returnHisto;
531 }