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) :
  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();
  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() {
 118   // initialize container obj
 119   FILE* fp0 = fopen(fFileName.Data(), "r");
 120   if (fp0 == NULL) {
 121     return;
 122   }
 123   else {
 124     fclose(fp0);
 125   }
 126   fPhysObj = new AliITSOnlineSPDphys(fFileName.Data());
 127 }
 128
 129 void AliITSOnlineSPDphysAnalyzer::SetParam(const Char_t *pname, const Char_t *pval) {
 130   // set a parameter
 131   TString name = pname;
 132   TString val = pval;
 133   //  printf("Setting Param %s  to %s\n",name.Data(),val.Data());
 134   if (name.CompareTo("MistakeProbabilityNoisy")==0) {
 135     Double_t mistakeProbabilityNoisy = val.Atof();
 136     fThreshNoisy = mistakeProbabilityNoisy/(20*6*10*32*256);
 137   }
 138   else if (name.CompareTo("MistakeProbabilityDead")==0) {
 139     Double_t mistakeProbabilityDead = val.Atof();
 140     fThreshDead = mistakeProbabilityDead/(20*6*10*32*256);
 141   }
 142   else if (name.CompareTo("fMinEventsForNoisy")==0) {
 143     fMinEventsForNoisy = val.Atoi();
 144   }
 145   else if (name.CompareTo("fMinEventsForDead")==0) {
 146     fMinEventsForDead = val.Atoi();
 147   }
 148   else if (name.CompareTo("fDefinitelyNoisyRatio")==0) {
 149     fDefinitelyNoisyRatio = val.Atof();
 150   }
 151   else if (name.CompareTo("fMinNrEqHitsForDeadChips")==0) {
 152     fMinNrEqHitsForDeadChips = val.Atof();
 153   }
 154   else if (name.CompareTo("fRatioToMeanForInefficientChip")==0) {
 155     fRatioToMeanForInefficientChip = val.Atof();
 156   }
 157   else {
 158     Error("AliITSOnlineSPDphysAnalyzer::SetParam","Parameter %s in configuration file unknown.",name.Data());
 159   }
 160 }
 161
 162 void AliITSOnlineSPDphysAnalyzer::ReadParamsFromLocation(const Char_t *dirName) {
 163   // opens file (default name) in dir dirName and reads parameters from it
 164   TString paramsFileName = Form("%s/physics_params.txt",dirName);
 165   ifstream paramsFile;
 166   paramsFile.open(paramsFileName, ifstream::in);
 167   if (paramsFile.fail()) {
 168     printf("No config file (%s) present. Using default tuning parameters.\n",paramsFileName.Data());
 169   }
 170   else {
 171     while(1) {
 172       Char_t paramN[50];
 173       Char_t paramV[50];
 174       paramsFile >> paramN;
 175       if (paramsFile.eof()) break;
 176       paramsFile >> paramV;
 177       SetParam(paramN,paramV);
 178       if (paramsFile.eof()) break;
 179     }
 180     paramsFile.close();
 181   }
 182 }
 183
 184
 185 UInt_t AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels() {
 186   // process noisy pixel data , returns number of chips with enough statistics
 187   if (fPhysObj==NULL) {
 188     Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","No data!");
 189     return 0;
 190   }
 191   // do we have enough events to even try the algorithm?
 192   if (GetNrEvents() < fMinEventsForNoisy) {
 193     Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Nr events (%d) < fMinEventsForNoisy (%d)!",GetNrEvents(),fMinEventsForNoisy);
 194     return 0;
 195   }
 196   // handler should be initialized
 197   if (fHandler==NULL) {
 198     Error("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Calibration handler is not initialized!");
 199     return 0;
 200   }
 201
 202   UInt_t nrEnoughStatChips = 0;
 203
 204   for (UInt_t hs=0; hs<6; hs++) {
 205     for (UInt_t chip=0; chip<10; chip++) {
 206
 207       UInt_t nrPixels = 0;
 208       UInt_t nrChipHits = 0;
 209       UInt_t nrMostHits = 0;
 210       for (UInt_t col=0; col<32; col++) {
 211         for (UInt_t row=0; row<256; row++) {
 212           UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
 213           nrChipHits += nrHits;
 214           //      if (nrHits>0) nrPixels++; // don't include pixels that might be dead
 215           nrPixels++;
 216           if (nrHits>fDefinitelyNoisyRatio*GetNrEvents()) {
 217             fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
 218             nrPixels--;
 219             nrChipHits-=nrHits;
 220           }
 221           else {
 222             if (nrMostHits<nrHits) nrMostHits=nrHits;
 223           }
 224         }
 225       }
 226
 227       if (nrChipHits>0) { // otherwise there are for sure no noisy
 228         // Binomial with n events and probability p for pixel hit
 229         UInt_t n = GetNrEvents();
 230         if (nrPixels>0 && n>0) {
 231
 232           Double_t p = (Double_t)nrChipHits/nrPixels/n;
 233
 234           // Bin(n,k=0):
 235           Double_t bin = pow((Double_t)(1-p),(Double_t)n);
 236           // Bin(n,k)
 237           UInt_t k=1;
 238           while ((bin>fThreshNoisy || k<n*p) && k<=n) {
 239             k++;
 240             bin = bin*(n-k+1)/k*p/(1-p);
 241           }
 242
 243           // can we find noisy pixels...?
 244           if (k<=n) {
 245             //      printf("eq %d , hs %d , chip %d : Noisy level = %d\n",GetEqNr(),hs,chip,k);
 246             nrEnoughStatChips++;
 247             // add noisy pixels to handler
 248             UInt_t noiseLimit=k;
 249             if (nrMostHits>=noiseLimit) {
 250               for (UInt_t col=0; col<32; col++) {
 251                 for (UInt_t row=0; row<256; row++) {
 252                   UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
 253                   if (nrHits >= noiseLimit) {
 254                     fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
 255                   }
 256                 }
 257               }
 258             }
 259           }
 260         }
 261
 262       }
 263
 264     } // for chip
 265   } // for hs
 266
 267   return nrEnoughStatChips;
 268 }
 269
 270
 271 UInt_t AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels() {
 272   // process dead pixel data , returns number of chips with enough statistics
 273   if (fPhysObj==NULL) {
 274     Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","No data!");
 275     return 0;
 276   }
 277
 278   // do we have enough events to even try the algorithm?
 279   if (GetNrEvents() < fMinEventsForDead) {
 280     Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Nr events (%d) < fMinEventsForDead (%d)!",GetNrEvents(),fMinEventsForDead);
 281     return 0;
 282   }
 283   // handler should be initialized
 284   if (fHandler==NULL) {
 285     Error("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Calibration handler is not initialized!");
 286     return 0;
 287   }
 288
 289   AliITSIntMap* possiblyDead  = new AliITSIntMap();
 290   AliITSIntMap* possiblyIneff = new AliITSIntMap();
 291
 292   fNrEnoughStatChips = 0;
 293   fNrDeadChips = 0;
 294   fNrInefficientChips = 0;
 295   UInt_t nrPossiblyDeadChips = 0;
 296   fNrEqHits = 0;
 297
 298
 299   for (UInt_t hs=0; hs<6; hs++) {
 300     if (!fHandler->IsActiveHS(GetEqNr(),hs)) {
 301       fNrDeadChips+=10;
 302     }
 303     else {
 304       for (UInt_t chip=0; chip<10; chip++) {
 305         if (!fHandler->IsActiveChip(GetEqNr(),hs,chip)) {
 306           fNrDeadChips++;
 307         }
 308         else {
 309           // perform search for individual dead pixels...
 310           Bool_t good=kFALSE;
 311
 312           UInt_t nrPossiblyDeadPixels = 0;
 313           UInt_t nrPixels = 0;
 314           UInt_t nrChipHits = 0;
 315           for (UInt_t col=0; col<32; col++) {
 316             for (UInt_t row=0; row<256; row++) {
 317               UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
 318               nrChipHits += nrHits;
 319               if (!fHandler->IsPixelNoisy(GetEqNr(),hs,chip,col,row)) {
 320                 // don't include noisy pixels
 321                 nrPixels++;
 322                 if (nrHits==0) {
 323                   nrPossiblyDeadPixels++;
 324                 }
 325                 else {
 326                   fHandler->UnSetDeadPixel(GetEqNr(),hs,chip,col,row); // unset (no action unless dead before)
 327                 }
 328               }
 329               else {
 330                 nrChipHits -= nrHits; // this is needed when running offline (online nrHits should be 0 already)
 331               }
 332             }
 333           }
 334           fNrEqHits+=nrChipHits;
 335
 336           if (nrChipHits>0) {
 337             // make sure the chip is not flagged as dead
 338             fHandler->SetDeadChip(GetEqNr(),hs,chip,kFALSE);
 339           }
 340
 341           if (nrPossiblyDeadPixels==0) {
 342             // no need to see if we have enough statistics...
 343             fNrEnoughStatChips++;
 344             good=kTRUE;
 345             //  printf("%3d",good);
 346             //  if (chip==9) printf("\n");
 347             continue;
 348           }
 349
 350           if (nrChipHits==0) {
 351             nrPossiblyDeadChips++;
 352             possiblyDead->Insert(hs,chip);
 353             good=kFALSE;
 354             //  printf("%3d",good);
 355             //  if (chip==9) printf("\n");
 356             continue;
 357           }
 358
 359           // Binomial with n events and probability p for pixel hit
 360           UInt_t n = GetNrEvents();
 361           if (nrPixels>0 && n>0) {
 362
 363             Double_t p = (Double_t)nrChipHits/nrPixels/n;
 364
 365             // probability of falsely assigning a dead pixel
 366             Double_t falselyDeadProb = pow((Double_t)(1-p),(Double_t)n);
 367             //      printf("falselyprob=%e\n",falselyDeadProb);
 368
 369             // can we find dead pixels...?
 370             if (falselyDeadProb<fThreshDead) {
 371               fNrEnoughStatChips++;
 372               good=kTRUE;
 373               // add dead pixels to handler
 374               for (UInt_t col=0; col<32; col++) {
 375                 for (UInt_t row=0; row<256; row++) {
 376                   UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
 377                   if (nrHits==0) {
 378                     if (!fHandler->IsPixelNoisy(GetEqNr(),hs,chip,col,row)) {
 379                       // don't include noisy pixels
 380                       fHandler->SetDeadPixel(GetEqNr(),hs,chip,col,row);
 381                     }
 382                   }
 383                 }
 384               }
 385             }
 386             if (!good) {
 387               // this might be an inefficient chip
 388               possiblyIneff->Insert(hs*10+chip,nrChipHits);
 389             }
 390
 391           }
 392           else {
 393             if (n>0) {
 394               // this is a completely noisy chip... put in category enough stat
 395               fNrEnoughStatChips++;
 396               good=kTRUE;
 397             }
 398           }
 399
 400           //      printf("%3d",good);
 401           //      if (chip==9) printf("\n");
 402
 403         }
 404       } // for chip
 405     }
 406   } // for hs
 407
 408
 409   Int_t key,val;
 410
 411   // dead chips?
 412   if (fNrEqHits>fMinNrEqHitsForDeadChips) {
 413     while (possiblyDead->Pop(key,val)) {
 414       fHandler->SetDeadChip(GetEqNr(),key,val,kFALSE);
 415     }
 416     fNrDeadChips+=nrPossiblyDeadChips;
 417   }
 418   delete possiblyDead;
 419
 420   // inefficient chips?
 421   while (possiblyIneff->Pop(key,val)) {
 422     if (val<fNrEqHits/60*fRatioToMeanForInefficientChip) {
 423       fNrInefficientChips++;
 424     }
 425   }
 426   delete possiblyIneff;
 427
 428
 429   fbDeadProcessed = kTRUE;
 430
 431   return fNrEnoughStatChips;
 432 }
 433
 434
 435 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEnoughStatChips() {
 436   // returns nr of enough stat chips
 437   if (!fbDeadProcessed) ProcessDeadPixels();
 438   return fNrEnoughStatChips;
 439 }
 440 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrDeadChips() {
 441   // returns nr of dead chips
 442   if (!fbDeadProcessed) ProcessDeadPixels();
 443   return fNrDeadChips;
 444 }
 445 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrInefficientChips() {
 446   // returns nr of inefficient chips
 447   if (!fbDeadProcessed) ProcessDeadPixels();
 448   return fNrInefficientChips;
 449 }
 450 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrNeedsMoreStatChips() {
 451   // returns nr of needs more stat chips
 452   if (!fbDeadProcessed) ProcessDeadPixels();
 453   return 60-fNrEnoughStatChips-fNrDeadChips-fNrInefficientChips;
 454 }
 455
 456 UInt_t AliITSOnlineSPDphysAnalyzer::GetEqNr() const {
 457   // returns the eq nr of phys obj
 458   if (fPhysObj!=NULL) return fPhysObj->GetEqNr();
 459   else return 999;
 460 }
 461
 462 UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEvents() const {
 463   // returns the nr of events of phys obj
 464   if (fPhysObj!=NULL) return fPhysObj->GetNrEvents();
 465   else return 0;
 466 }
 467
 468 void AliITSOnlineSPDphysAnalyzer::Exponent(Double_t &val, Int_t &valExp) const {
 469   // put double in format with val and exp so that 1<val<10 - The actual value is val*10e(valExp)
 470   while (val>10) {
 471     val/=10;
 472     valExp++;
 473   }
 474   while (val<1) {
 475     val*=10;
 476     valExp--;
 477   }
 478 }
 479
 480
 481
 482 TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapTot() {
 483   // creates and returns a pointer to a hitmap histo (half sector style a la spdmood)
 484   if (fPhysObj==NULL) {
 485     Error("AliITSOnlineSPDphysAnalyzer::GetHitMapTot","No data!");
 486     return NULL;
 487   }
 488   TString histoname = Form("Eq %d",GetEqNr());
 489   TH2F* fHitMapTot = new TH2F(histoname.Data(),histoname.Data(),32*10,-0.5,32*10-0.5,256*6,-0.5,256*6-0.5);
 490   fHitMapTot->SetNdivisions(-10,"X");
 491   fHitMapTot->SetNdivisions(-006,"Y");
 492   fHitMapTot->SetTickLength(0,"X");
 493   fHitMapTot->SetTickLength(0,"Y");
 494   fHitMapTot->GetXaxis()->SetLabelColor(gStyle->GetCanvasColor());
 495   fHitMapTot->GetYaxis()->SetLabelColor(gStyle->GetCanvasColor());
 496   for (UInt_t hs=0; hs<6; hs++) {
 497     for (UInt_t chipNr=0; chipNr<10; chipNr++) {
 498       for (UInt_t col=0; col<32; col++) {
 499         for (UInt_t row=0; row<256; row++) {
 500           fHitMapTot->Fill(chipNr*32+col,(5-hs)*256+row,fPhysObj->GetHits(hs,chipNr,col,row));
 501         }
 502       }
 503     }
 504   }
 505   return fHitMapTot;
 506 }
 507
 508 TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapChip(UInt_t hs, UInt_t chip) {
 509   // creates and returns a pointer to a hitmap histo (chip style a la spdmood)
 510   if (fPhysObj==NULL) {
 511     Error("AliITSOnlineSPDphysAnalyzer::GetHitMapChip","No data!");
 512     return NULL;
 513   }
 514
 515   TString histoName;
 516   TString histoTitle;
 517   histoName = Form("fChipHisto_%d_%d_%d", GetEqNr(), hs, chip);
 518   histoTitle = Form("Eq ID %d, Half Stave %d, Chip %d", GetEqNr(), hs, chip);
 519
 520   TH2F *returnHisto = new TH2F(histoName.Data(), histoTitle.Data(), 32, -0.5, 31.5, 256, -0.5, 255.5);
 521   returnHisto->SetMinimum(0);
 522   for (UInt_t col=0; col<32; col++) {
 523     for (UInt_t row=0; row<256; row++) {
 524       returnHisto->Fill(col,row,fPhysObj->GetHits(hs,chip,col,row));
 525     }
 526   }
 527
 528   return returnHisto;
 529 }