1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 ////////////////////////////////////////////////////////////
23 #include "AliITSOnlineSPDphysAnalyzer.h"
24 #include "AliITSOnlineSPDphys.h"
25 #include "AliITSOnlineCalibrationSPDhandler.h"
26 #include "AliITSRawStreamSPD.h"
27 #include "AliITSIntMap.h"
37 AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(const Char_t *fileName, AliITSOnlineCalibrationSPDhandler* handler) :
38 fFileName(fileName),fPhysObj(NULL),fHandler(handler),
39 fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
40 fNrEqHits(0),fbDeadProcessed(kFALSE),
41 fThreshNoisy(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
42 fMinEventsForNoisy(10000),fMinEventsForDead(10000),
43 fDefinitelyNoisyRatio(0.1),
44 fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)//!!!
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(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
55 fMinEventsForNoisy(10000),fMinEventsForDead(10000),
56 fDefinitelyNoisyRatio(0.1),
57 fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)//!!!
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(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
68 fMinEventsForNoisy(10000),fMinEventsForDead(10000),
69 fDefinitelyNoisyRatio(0.1),
70 fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)//!!!
72 // copy constructor, only copies the filename and params (not the processed data)
73 fFileName=handle.fFileName;
74 fThreshNoisy = handle.fThreshNoisy;
75 fThreshNoisyExp = handle.fThreshNoisyExp;
76 fThreshDead = handle.fThreshDead;
77 fThreshDeadExp = handle.fThreshDeadExp;
78 fMinEventsForNoisy = handle.fMinEventsForNoisy;
79 fMinEventsForDead = handle.fMinEventsForDead;
80 fDefinitelyNoisyRatio = handle.fDefinitelyNoisyRatio;
81 fMinNrEqHitsForDeadChips = handle.fMinNrEqHitsForDeadChips;
82 fRatioToMeanForInefficientChip = handle.fRatioToMeanForInefficientChip;
87 AliITSOnlineSPDphysAnalyzer::~AliITSOnlineSPDphysAnalyzer() {
89 if (fPhysObj!=NULL) delete fPhysObj;
92 AliITSOnlineSPDphysAnalyzer& AliITSOnlineSPDphysAnalyzer::operator=(const AliITSOnlineSPDphysAnalyzer& handle) {
93 // assignment operator, only copies the filename and params (not the processed data)
95 if (fPhysObj!=NULL) delete fPhysObj;
97 fFileName=handle.fFileName;
98 fThreshNoisy = handle.fThreshNoisy;
99 fThreshNoisyExp = handle.fThreshNoisyExp;
100 fThreshDead = handle.fThreshDead;
101 fThreshDeadExp = handle.fThreshDeadExp;
102 fMinEventsForNoisy = handle.fMinEventsForNoisy;
103 fMinEventsForDead = handle.fMinEventsForDead;
104 fDefinitelyNoisyRatio = handle.fDefinitelyNoisyRatio;
105 fMinNrEqHitsForDeadChips = handle.fMinNrEqHitsForDeadChips;
106 fRatioToMeanForInefficientChip = handle.fRatioToMeanForInefficientChip;
110 fNrEnoughStatChips = 0;
112 fNrInefficientChips = 0;
114 fbDeadProcessed = kFALSE;
121 void AliITSOnlineSPDphysAnalyzer::Init() {
122 // initialize container obj
123 FILE* fp0 = fopen(fFileName.Data(), "r");
130 fPhysObj = new AliITSOnlineSPDphys(fFileName.Data());
133 void AliITSOnlineSPDphysAnalyzer::SetParam(const Char_t *pname, const Char_t *pval) {
135 TString name = pname;
137 // printf("Setting Param %s to %s\n",name.Data(),val.Data());
138 if (name.CompareTo("MistakeProbabilityNoisy")==0) {
139 Double_t mistakeProbabilityNoisy = val.Atof();
140 fThreshNoisy = mistakeProbabilityNoisy;
142 Exponent(fThreshNoisy,fThreshNoisyExp);
143 fThreshNoisy/=(20*6*10*32*256);
144 Exponent(fThreshNoisy,fThreshNoisyExp);
146 else if (name.CompareTo("MistakeProbabilityDead")==0) {
147 Double_t mistakeProbabilityDead = val.Atof();
148 fThreshDead = mistakeProbabilityDead;
150 Exponent(fThreshDead,fThreshDeadExp);
151 fThreshDead/=(20*6*10*32*256);
152 Exponent(fThreshDead,fThreshDeadExp);
154 else if (name.CompareTo("fMinEventsForNoisy")==0) {
155 fMinEventsForNoisy = val.Atoi();
157 else if (name.CompareTo("fMinEventsForDead")==0) {
158 fMinEventsForDead = val.Atoi();
160 else if (name.CompareTo("fDefinitelyNoisyRatio")==0) {
161 fDefinitelyNoisyRatio = val.Atof();
163 else if (name.CompareTo("fMinNrEqHitsForDeadChips")==0) {
164 fMinNrEqHitsForDeadChips = val.Atof();
166 else if (name.CompareTo("fRatioToMeanForInefficientChip")==0) {
167 fRatioToMeanForInefficientChip = val.Atof();
170 Error("AliITSOnlineSPDphysAnalyzer::SetParam","Parameter %s in configuration file unknown.",name.Data());
174 void AliITSOnlineSPDphysAnalyzer::ReadParamsFromLocation(const Char_t *dirName) {
175 // opens file (default name) in dir dirName and reads parameters from it
176 TString paramsFileName = Form("%s/physics_params.txt",dirName);
178 paramsFile.open(paramsFileName, ifstream::in);
179 if (paramsFile.fail()) {
180 printf("No config file (%s) present. Using default tuning parameters.\n",paramsFileName.Data());
186 paramsFile >> paramN;
187 if (paramsFile.eof()) break;
188 paramsFile >> paramV;
189 SetParam(paramN,paramV);
190 if (paramsFile.eof()) break;
197 UInt_t AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels() {
198 // process noisy pixel data , returns number of chips with enough statistics
199 if (fPhysObj==NULL) {
200 Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","No data!");
203 // do we have enough events to even try the algorithm?
204 if (GetNrEvents() < fMinEventsForNoisy) {
205 Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Nr events (%d) < fMinEventsForNoisy (%d)!",GetNrEvents(),fMinEventsForNoisy);
208 // handler should be initialized
209 if (fHandler==NULL) {
210 Error("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Calibration handler is not initialized!");
214 UInt_t nrEnoughStatChips = 0;
216 for (UInt_t hs=0; hs<6; hs++) {
217 for (UInt_t chip=0; chip<10; chip++) {
220 UInt_t nrChipHits = 0;
221 UInt_t nrMostHits = 0;
222 for (UInt_t col=0; col<32; col++) {
223 for (UInt_t row=0; row<256; row++) {
224 UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
225 nrChipHits += nrHits;
226 // if (nrHits>0) nrPixels++; // don't include pixels that might be dead
228 if (nrHits>fDefinitelyNoisyRatio*GetNrEvents()) {
229 fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
234 if (nrMostHits<nrHits) nrMostHits=nrHits;
239 if (nrChipHits>0) { // otherwise there are for sure no noisy
240 // Binomial with n events and probability p for pixel hit
241 UInt_t n = GetNrEvents();
242 if (nrPixels>0 && n>0) {
244 Double_t p = (Double_t)nrChipHits/nrPixels/n;
249 for (UInt_t i=0; i<n; i++) {
251 Exponent(bin,binExp);
255 while (((binExp>fThreshNoisyExp || (binExp==fThreshNoisyExp && bin>fThreshNoisy)) || k<n*p) && k<=n) {
257 bin = bin*(n-k+1)/k*p/(1-p);
258 Exponent(bin,binExp);
261 // can we find noisy pixels...?
263 // printf("eq %d , hs %d , chip %d : Noisy level = %d\n",GetEqNr(),hs,chip,k);
265 // add noisy pixels to handler
267 if (nrMostHits>=noiseLimit) {
268 for (UInt_t col=0; col<32; col++) {
269 for (UInt_t row=0; row<256; row++) {
270 UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
271 if (nrHits >= noiseLimit) {
272 fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
285 return nrEnoughStatChips;
289 UInt_t AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels() {
290 // process dead pixel data , returns number of chips with enough statistics
291 if (fPhysObj==NULL) {
292 Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","No data!");
295 // do we have enough events to even try the algorithm?
296 if (GetNrEvents() < fMinEventsForDead) {
297 Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Nr events (%d) < fMinEventsForDead (%d)!",GetNrEvents(),fMinEventsForDead);
300 // handler should be initialized
301 if (fHandler==NULL) {
302 Error("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Calibration handler is not initialized!");
306 AliITSOnlineCalibrationSPDhandler *deadChipHandler = new AliITSOnlineCalibrationSPDhandler();
307 AliITSIntMap* possiblyIneff = new AliITSIntMap();
309 fNrEnoughStatChips = 0;
311 fNrInefficientChips = 0;
312 UInt_t nrPossiblyDeadChips = 0;
315 for (UInt_t hs=0; hs<6; hs++) {
316 for (UInt_t chip=0; chip<10; chip++) {
319 AliITSOnlineCalibrationSPDhandler *deadPixelHandler = new AliITSOnlineCalibrationSPDhandler();
321 UInt_t nrPossiblyDeadPixels = 0;
323 UInt_t nrChipHits = 0;
324 for (UInt_t col=0; col<32; col++) {
325 for (UInt_t row=0; row<256; row++) {
326 UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
327 nrChipHits += nrHits;
328 if (!fHandler->IsPixelNoisy(GetEqNr(),hs,chip,col,row)) {
329 // don't include noisy pixels
332 nrPossiblyDeadPixels++;
333 deadPixelHandler->SetDeadPixel(GetEqNr(),hs,chip,col,row);
337 nrChipHits -= nrHits; // this is needed when running offline (online nrHits should be 0 already)
341 fNrEqHits+=nrChipHits;
343 // check all pixels that were declared dead from before...
344 UInt_t nrDeadBefore = fHandler->GetNrDeadC(GetEqNr(),hs,chip);
345 UInt_t nrNoisyNow = fHandler->GetNrNoisyC(GetEqNr(),hs,chip);
346 if (nrDeadBefore+nrNoisyNow==8192) {
348 fHandler->UnSetDeadChip(GetEqNr(),hs,chip);
352 AliITSOnlineCalibrationSPDhandler *tmpHand = new AliITSOnlineCalibrationSPDhandler();
353 for (UInt_t index=0; index<nrDeadBefore; index++) {
354 UInt_t col = fHandler->GetDeadColAtC(GetEqNr(),hs,chip,index);
355 UInt_t row = fHandler->GetDeadRowAtC(GetEqNr(),hs,chip,index);
356 if (fPhysObj->GetHits(hs,chip,col,row)>0) {
357 // fHandler->UnSetDeadPixel(GetEqNr(),hs,chip,col,row);
358 tmpHand->SetDeadPixel(GetEqNr(),hs,chip,col,row);
361 UInt_t nrToRemove = tmpHand->GetNrDead();
362 for (UInt_t index=0; index<nrToRemove; index++) {
363 fHandler->UnSetDeadPixel(GetEqNr(),hs,chip,tmpHand->GetDeadColAtC(GetEqNr(),hs,chip,index),tmpHand->GetDeadRowAtC(GetEqNr(),hs,chip,index));
371 if (nrPossiblyDeadPixels==0) {
372 // no need to see if we have enough statistics...
373 fNrEnoughStatChips++;
375 delete deadPixelHandler;
378 if (chip==9) printf("\n");
383 nrPossiblyDeadChips++;
384 deadChipHandler->SetDeadChip(GetEqNr(),hs,chip);
385 delete deadPixelHandler;
388 if (chip==9) printf("\n");
392 // Binomial with n events and probability p for pixel hit
393 UInt_t n = GetNrEvents();
394 if (nrPixels>0 && n>0) {
396 Double_t p = (Double_t)nrChipHits/nrPixels/n;
398 // probability of falsely assigning a dead pixel
399 Double_t falselyDeadProb = 1;
400 Int_t falselyDeadProbExp = 0;
401 for (UInt_t i=0; i<n; i++) {
402 falselyDeadProb*=(1-p);
403 Exponent(falselyDeadProb,falselyDeadProbExp);
406 // can we find dead pixels...?
407 if (falselyDeadProbExp<fThreshDeadExp || (falselyDeadProbExp==fThreshDeadExp && falselyDeadProb<fThreshDead)) {
408 fNrEnoughStatChips++;
410 // add dead pixels to handler
411 fHandler->AddDeadFrom(deadPixelHandler);
414 // this might be an inefficient chip
415 possiblyIneff->Insert(hs*10+chip,nrChipHits);
421 // this is a completely noisy chip... put in category enough stat
422 fNrEnoughStatChips++;
427 delete deadPixelHandler;
430 if (chip==9) printf("\n");
437 if (fNrEqHits>fMinNrEqHitsForDeadChips) {
438 fHandler->AddDeadFrom(deadChipHandler);
439 fNrDeadChips=nrPossiblyDeadChips;
444 delete deadChipHandler;
447 // inefficient chips?
449 while (possiblyIneff->Pop(key,val)) {
450 if (val<fNrEqHits/60*fRatioToMeanForInefficientChip) {
451 fNrInefficientChips++;
454 delete possiblyIneff;
457 fbDeadProcessed = kTRUE;
459 return fNrEnoughStatChips;
463 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEnoughStatChips() {
464 // returns nr of enough stat chips
465 if (!fbDeadProcessed) ProcessDeadPixels();
466 return fNrEnoughStatChips;
468 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrDeadChips() {
469 // returns nr of dead chips
470 if (!fbDeadProcessed) ProcessDeadPixels();
473 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrInefficientChips() {
474 // returns nr of inefficient chips
475 if (!fbDeadProcessed) ProcessDeadPixels();
476 return fNrInefficientChips;
478 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrNeedsMoreStatChips() {
479 // returns nr of needs more stat chips
480 if (!fbDeadProcessed) ProcessDeadPixels();
481 return 60-fNrEnoughStatChips-fNrDeadChips-fNrInefficientChips;
484 UInt_t AliITSOnlineSPDphysAnalyzer::GetEqNr() const {
485 // returns the eq nr of phys obj
486 if (fPhysObj!=NULL) return fPhysObj->GetEqNr();
490 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEvents() const {
491 // returns the nr of events of phys obj
492 if (fPhysObj!=NULL) return fPhysObj->GetNrEvents();
496 void AliITSOnlineSPDphysAnalyzer::Exponent(Double_t &val, Int_t &valExp) const {
497 // put double in format with val and exp so that 1<val<10 - The actual value is val*10e(valExp)
510 TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapTot() {
511 // creates and returns a pointer to a hitmap histo (half sector style a la spdmood)
512 if (fPhysObj==NULL) {
513 Error("AliITSOnlineSPDphysAnalyzer::GetHitMapTot","No data!");
516 TString histoname = Form("Eq %d",GetEqNr());
517 TH2F* fHitMapTot = new TH2F(histoname.Data(),histoname.Data(),32*10,-0.5,32*10-0.5,256*6,-0.5,256*6-0.5);
518 fHitMapTot->SetNdivisions(-10,"X");
519 fHitMapTot->SetNdivisions(-006,"Y");
520 fHitMapTot->SetTickLength(0,"X");
521 fHitMapTot->SetTickLength(0,"Y");
522 fHitMapTot->GetXaxis()->SetLabelColor(gStyle->GetCanvasColor());
523 fHitMapTot->GetYaxis()->SetLabelColor(gStyle->GetCanvasColor());
524 for (UInt_t hs=0; hs<6; hs++) {
525 for (UInt_t chipNr=0; chipNr<10; chipNr++) {
526 for (UInt_t col=0; col<32; col++) {
527 for (UInt_t row=0; row<256; row++) {
528 fHitMapTot->Fill(chipNr*32+col,(5-hs)*256+row,fPhysObj->GetHits(hs,chipNr,col,row));
536 TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapChip(UInt_t hs, UInt_t chip) {
537 // creates and returns a pointer to a hitmap histo (chip style a la spdmood)
538 if (fPhysObj==NULL) {
539 Error("AliITSOnlineSPDphysAnalyzer::GetHitMapChip","No data!");
545 histoName = Form("fChipHisto_%d_%d_%d", GetEqNr(), hs, chip);
546 histoTitle = Form("Eq ID %d, Half Stave %d, Chip %d", GetEqNr(), hs, chip);
548 TH2F *returnHisto = new TH2F(histoName.Data(), histoTitle.Data(), 32, -0.5, 31.5, 256, -0.5, 255.5);
549 returnHisto->SetMinimum(0);
550 for (UInt_t col=0; col<32; col++) {
551 for (UInt_t row=0; row<256; row++) {
552 returnHisto->Fill(col,row,fPhysObj->GetHits(hs,chip,col,row));