[u/mrichter/AliRoot.git] / ITS / AliITSOnlineSPDphysAnalyzer.cxx

/**************************************************************************
 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

////////////////////////////////////////////////////////////
// Author: Henrik Tydesjo                                 //
// This class is used in the detector algorithm framework //
// to process the data stored in special container files  //
// (see AliITSOnlineSPDphys).                             //
////////////////////////////////////////////////////////////

#include "AliITSOnlineSPDphysAnalyzer.h"
#include "AliITSOnlineSPDphys.h"
#include "AliITSOnlineCalibrationSPDhandler.h"
#include "AliITSRawStreamSPD.h"
#include "AliITSIntMap.h"
#include <TStyle.h>
#include <TMath.h>
#include <TF1.h>
#include <TGraph.h>
#include <TH2F.h>
#include <TError.h>
#include <iostream>
#include <fstream>

AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(const Char_t *fileName, AliITSOnlineCalibrationSPDhandler* handler) :
  fFileName(fileName),fPhysObj(NULL),fHandler(handler),
  fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
  fNrEqHits(0),fbDeadProcessed(kFALSE),
  fThreshNoisy(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
  fMinEventsForNoisy(10000),fMinEventsForDead(10000),
  fDefinitelyNoisyRatio(0.1),
  fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)//!!!
{
  // constructor  
  Init();
}

AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(AliITSOnlineSPDphys* physObj, AliITSOnlineCalibrationSPDhandler* handler) :
  fFileName("test.root"),fPhysObj(NULL),fHandler(handler),
  fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
  fNrEqHits(0),fbDeadProcessed(kFALSE),
  fThreshNoisy(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
  fMinEventsForNoisy(10000),fMinEventsForDead(10000),
  fDefinitelyNoisyRatio(0.1),
  fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)//!!!
{
  // alt constructor
  fPhysObj = physObj;
}

AliITSOnlineSPDphysAnalyzer::AliITSOnlineSPDphysAnalyzer(const AliITSOnlineSPDphysAnalyzer& handle) :
  fFileName("test.root"),fPhysObj(NULL),fHandler(NULL),
  fNrEnoughStatChips(0),fNrDeadChips(0),fNrInefficientChips(0),
  fNrEqHits(0),fbDeadProcessed(kFALSE),
  fThreshNoisy(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
  fMinEventsForNoisy(10000),fMinEventsForDead(10000),
  fDefinitelyNoisyRatio(0.1),
  fMinNrEqHitsForDeadChips(60000),fRatioToMeanForInefficientChip(0.1)//!!!
{
  // copy constructor, only copies the filename and params (not the processed data)
  fFileName=handle.fFileName;
  fThreshNoisy = handle.fThreshNoisy;
  fThreshNoisyExp = handle.fThreshNoisyExp;
  fThreshDead = handle.fThreshDead;
  fThreshDeadExp = handle.fThreshDeadExp;
  fMinEventsForNoisy = handle.fMinEventsForNoisy;
  fMinEventsForDead = handle.fMinEventsForDead;
  fDefinitelyNoisyRatio = handle.fDefinitelyNoisyRatio;
  fMinNrEqHitsForDeadChips = handle.fMinNrEqHitsForDeadChips;
  fRatioToMeanForInefficientChip = handle.fRatioToMeanForInefficientChip;

  Init();
}

AliITSOnlineSPDphysAnalyzer::~AliITSOnlineSPDphysAnalyzer() {
  // destructor
  if (fPhysObj!=NULL) delete fPhysObj;
}

AliITSOnlineSPDphysAnalyzer& AliITSOnlineSPDphysAnalyzer::operator=(const AliITSOnlineSPDphysAnalyzer& handle) {
  // assignment operator, only copies the filename and params (not the processed data)
  if (this!=&handle) {
    if (fPhysObj!=NULL) delete fPhysObj;
    
    fFileName=handle.fFileName;
    fThreshNoisy = handle.fThreshNoisy;
    fThreshNoisyExp = handle.fThreshNoisyExp;
    fThreshDead = handle.fThreshDead;
    fThreshDeadExp = handle.fThreshDeadExp;
    fMinEventsForNoisy = handle.fMinEventsForNoisy;
    fMinEventsForDead = handle.fMinEventsForDead;
    fDefinitelyNoisyRatio = handle.fDefinitelyNoisyRatio;
    fMinNrEqHitsForDeadChips = handle.fMinNrEqHitsForDeadChips;
    fRatioToMeanForInefficientChip = handle.fRatioToMeanForInefficientChip;

    fPhysObj = NULL;
    fHandler = NULL;
    fNrEnoughStatChips = 0;
    fNrDeadChips = 0;
    fNrInefficientChips = 0;
    fNrEqHits = 0;
    fbDeadProcessed = kFALSE;

    Init();    
  }
  return *this;
}

void AliITSOnlineSPDphysAnalyzer::Init() {
  // initialize container obj
  FILE* fp0 = fopen(fFileName.Data(), "r");
  if (fp0 == NULL) {
    return;
  }
  else {
    fclose(fp0);
  }
  fPhysObj = new AliITSOnlineSPDphys(fFileName.Data());
}

void AliITSOnlineSPDphysAnalyzer::SetParam(const Char_t *pname, const Char_t *pval) {
  // set a parameter
  TString name = pname;
  TString val = pval;
  //  printf("Setting Param %s  to %s\n",name.Data(),val.Data());
  if (name.CompareTo("MistakeProbabilityNoisy")==0) {
    Double_t mistakeProbabilityNoisy = val.Atof();
    fThreshNoisy = mistakeProbabilityNoisy;
    fThreshNoisyExp = 0;
    Exponent(fThreshNoisy,fThreshNoisyExp);
    fThreshNoisy/=(20*6*10*32*256);
    Exponent(fThreshNoisy,fThreshNoisyExp);
  }
  else if (name.CompareTo("MistakeProbabilityDead")==0) {
    Double_t mistakeProbabilityDead = val.Atof();
    fThreshDead = mistakeProbabilityDead;
    fThreshDeadExp = 0;
    Exponent(fThreshDead,fThreshDeadExp);
    fThreshDead/=(20*6*10*32*256);
    Exponent(fThreshDead,fThreshDeadExp);
  }
  else if (name.CompareTo("fMinEventsForNoisy")==0) {
    fMinEventsForNoisy = val.Atoi();
  }
  else if (name.CompareTo("fMinEventsForDead")==0) {
    fMinEventsForDead = val.Atoi();
  }
  else if (name.CompareTo("fDefinitelyNoisyRatio")==0) {
    fDefinitelyNoisyRatio = val.Atof();
  }
  else if (name.CompareTo("fMinNrEqHitsForDeadChips")==0) {
    fMinNrEqHitsForDeadChips = val.Atof();
  }
  else if (name.CompareTo("fRatioToMeanForInefficientChip")==0) {
    fRatioToMeanForInefficientChip = val.Atof();
  }
  else {
    Error("AliITSOnlineSPDphysAnalyzer::SetParam","Parameter %s in configuration file unknown.",name.Data());
  }
}

void AliITSOnlineSPDphysAnalyzer::ReadParamsFromLocation(const Char_t *dirName) {
  // opens file (default name) in dir dirName and reads parameters from it
  TString paramsFileName = Form("%s/physics_params.txt",dirName);
  ifstream paramsFile;
  paramsFile.open(paramsFileName, ifstream::in);
  if (paramsFile.fail()) {
    printf("No config file (%s) present. Using default tuning parameters.\n",paramsFileName.Data());
  }
  else {
    while(1) {
      Char_t paramN[50];
      Char_t paramV[50];
      paramsFile >> paramN;
      if (paramsFile.eof()) break;
      paramsFile >> paramV;
      SetParam(paramN,paramV);
      if (paramsFile.eof()) break;
    }
    paramsFile.close();
  }
}


UInt_t AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels() {
  // process noisy pixel data , returns number of chips with enough statistics
  if (fPhysObj==NULL) {
    Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","No data!");
    return 0;
  }
  // do we have enough events to even try the algorithm?
  if (GetNrEvents() < fMinEventsForNoisy) {
    Warning("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Nr events (%d) < fMinEventsForNoisy (%d)!",GetNrEvents(),fMinEventsForNoisy);
    return 0;
  }
  // handler should be initialized
  if (fHandler==NULL) {
    Error("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Calibration handler is not initialized!");
    return 0;
  }
  
  UInt_t nrEnoughStatChips = 0;

  for (UInt_t hs=0; hs<6; hs++) {
    for (UInt_t chip=0; chip<10; chip++) {

      UInt_t nrPixels = 0;
      UInt_t nrChipHits = 0;
      UInt_t nrMostHits = 0;
      for (UInt_t col=0; col<32; col++) {
	for (UInt_t row=0; row<256; row++) {
	  UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
	  nrChipHits += nrHits;
	  //	  if (nrHits>0) nrPixels++; // don't include pixels that might be dead
	  nrPixels++;
	  if (nrHits>fDefinitelyNoisyRatio*GetNrEvents()) {
	    fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
	    nrPixels--;
	    nrChipHits-=nrHits;
	  }
	  else {
	    if (nrMostHits<nrHits) nrMostHits=nrHits;
	  }
	}
      }

      if (nrChipHits>0) { // otherwise there are for sure no noisy
	// Binomial with n events and probability p for pixel hit
	UInt_t n = GetNrEvents();
	if (nrPixels>0 && n>0) {

	  Double_t p = (Double_t)nrChipHits/nrPixels/n;

	  Double_t bin = 1;
	  Int_t binExp = 0;
	  // Bin(n,k=0):
	  for (UInt_t i=0; i<n; i++) {
	    bin*=(1-p);
	    Exponent(bin,binExp);
	  }
	  // Bin(n,k)
	  UInt_t k=1;
	  while (((binExp>fThreshNoisyExp || (binExp==fThreshNoisyExp && bin>fThreshNoisy)) || k<n*p) && k<=n) {
	    k++;
	    bin = bin*(n-k+1)/k*p/(1-p);
	    Exponent(bin,binExp);
	  }
	  
	  // can we find noisy pixels...?
	  if (k<=n) {
	    //	    printf("eq %d , hs %d , chip %d : Noisy level = %d\n",GetEqNr(),hs,chip,k);
	    nrEnoughStatChips++;
	    // add noisy pixels to handler
	    UInt_t noiseLimit=k;
	    if (nrMostHits>=noiseLimit) {
	      for (UInt_t col=0; col<32; col++) {
		for (UInt_t row=0; row<256; row++) {
		  UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
		  if (nrHits >= noiseLimit) {
		    fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
		  }
		}
	      }
	    }
	  }
	}

      }

    } // for chip
  } // for hs

  return nrEnoughStatChips;
}


UInt_t AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels() {
  // process dead pixel data , returns number of chips with enough statistics
  if (fPhysObj==NULL) {
    Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","No data!");
    return 0;
  }
  // do we have enough events to even try the algorithm?
  if (GetNrEvents() < fMinEventsForDead) {
    Warning("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Nr events (%d) < fMinEventsForDead (%d)!",GetNrEvents(),fMinEventsForDead);
    return 0;
  }
  // handler should be initialized
  if (fHandler==NULL) {
    Error("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Calibration handler is not initialized!");
    return 0;
  }

  AliITSOnlineCalibrationSPDhandler *deadChipHandler = new AliITSOnlineCalibrationSPDhandler();
  AliITSIntMap* possiblyIneff = new AliITSIntMap();

  fNrEnoughStatChips = 0;
  fNrDeadChips = 0;
  fNrInefficientChips = 0;
  UInt_t nrPossiblyDeadChips = 0;
  fNrEqHits = 0;

  for (UInt_t hs=0; hs<6; hs++) {
    for (UInt_t chip=0; chip<10; chip++) {
      Bool_t good=kFALSE;

      AliITSOnlineCalibrationSPDhandler *deadPixelHandler = new AliITSOnlineCalibrationSPDhandler();

      UInt_t nrPossiblyDeadPixels = 0;
      UInt_t nrPixels = 0;
      UInt_t nrChipHits = 0;
      for (UInt_t col=0; col<32; col++) {
	for (UInt_t row=0; row<256; row++) {
	  UInt_t nrHits = fPhysObj->GetHits(hs,chip,col,row);
	  nrChipHits += nrHits;
	  if (!fHandler->IsPixelNoisy(GetEqNr(),hs,chip,col,row)) {
	    // don't include noisy pixels
	    nrPixels++;
	    if (nrHits==0) {
	      nrPossiblyDeadPixels++;
	      deadPixelHandler->SetDeadPixel(GetEqNr(),hs,chip,col,row);
	    }
	  }
	  else {
	    nrChipHits -= nrHits; // this is needed when running offline (online nrHits should be 0 already)
	  }
	}
      }
      fNrEqHits+=nrChipHits;

      // check all pixels that were declared dead from before...
      UInt_t nrDeadBefore = fHandler->GetNrDeadC(GetEqNr(),hs,chip);
      UInt_t nrNoisyNow   = fHandler->GetNrNoisyC(GetEqNr(),hs,chip);
      if (nrDeadBefore+nrNoisyNow==8192) {
	if (nrChipHits>0) {
	  fHandler->UnSetDeadChip(GetEqNr(),hs,chip);
	}
      }
      else {
	AliITSOnlineCalibrationSPDhandler *tmpHand = new AliITSOnlineCalibrationSPDhandler();
	for (UInt_t index=0; index<nrDeadBefore; index++) {
	  UInt_t col = fHandler->GetDeadColAtC(GetEqNr(),hs,chip,index);
	  UInt_t row = fHandler->GetDeadRowAtC(GetEqNr(),hs,chip,index);
	  if (fPhysObj->GetHits(hs,chip,col,row)>0) {
	    //	    fHandler->UnSetDeadPixel(GetEqNr(),hs,chip,col,row);
	    tmpHand->SetDeadPixel(GetEqNr(),hs,chip,col,row);
	  }
	}
	UInt_t nrToRemove = tmpHand->GetNrDead();
	for (UInt_t index=0; index<nrToRemove; index++) {
	  fHandler->UnSetDeadPixel(GetEqNr(),hs,chip,tmpHand->GetDeadColAtC(GetEqNr(),hs,chip,index),tmpHand->GetDeadRowAtC(GetEqNr(),hs,chip,index));
	}
	delete tmpHand;
      }


      if (nrPossiblyDeadPixels==0) {
	// no need to see if we have enough statistics...
	fNrEnoughStatChips++;
	good=kTRUE;
	delete deadPixelHandler;
	good=kTRUE;
	printf("%3d",good);
	if (chip==9) printf("\n");
	continue;
      }

      if (nrChipHits==0) {
	nrPossiblyDeadChips++;
	deadChipHandler->SetDeadChip(GetEqNr(),hs,chip);
	delete deadPixelHandler;
	good=kFALSE;
	printf("%3d",good);
	if (chip==9) printf("\n");
	continue;
      }

      // Binomial with n events and probability p for pixel hit
      UInt_t n = GetNrEvents();
      if (nrPixels>0 && n>0) {

	Double_t p = (Double_t)nrChipHits/nrPixels/n;

	// probability of falsely assigning a dead pixel
	Double_t falselyDeadProb = 1;
	Int_t falselyDeadProbExp = 0;
	for (UInt_t i=0; i<n; i++) {
	  falselyDeadProb*=(1-p);
	  Exponent(falselyDeadProb,falselyDeadProbExp);
	}

	// can we find dead pixels...?
	if (falselyDeadProbExp<fThreshDeadExp || (falselyDeadProbExp==fThreshDeadExp && falselyDeadProb<fThreshDead)) {
	  fNrEnoughStatChips++;
	  good=kTRUE;
	  // add dead pixels to handler
	  fHandler->AddDeadFrom(deadPixelHandler);
	}
	else {
	  // this might be an inefficient chip
	  possiblyIneff->Insert(hs*10+chip,nrChipHits);
	}

      }
      else {
	if (n>0) {
	  // this is a completely noisy chip... put in category enough stat
	  fNrEnoughStatChips++;
	  good=kTRUE;
	}
      }

      delete deadPixelHandler;

      printf("%3d",good);
      if (chip==9) printf("\n");

    } // for chip
  } // for hs


  // dead chips?
  if (fNrEqHits>fMinNrEqHitsForDeadChips) {
    fHandler->AddDeadFrom(deadChipHandler);
    fNrDeadChips=nrPossiblyDeadChips;
  }
  else {
    fNrDeadChips=0;
  }
  delete deadChipHandler;


  // inefficient chips?
  Int_t key,val;
  while (possiblyIneff->Pop(key,val)) {
    if (val<fNrEqHits/60*fRatioToMeanForInefficientChip) {
      fNrInefficientChips++;
    }
  }
  delete possiblyIneff;


  fbDeadProcessed = kTRUE;

  return fNrEnoughStatChips;
}


UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEnoughStatChips() {
  // returns nr of enough stat chips
  if (!fbDeadProcessed) ProcessDeadPixels();
  return fNrEnoughStatChips;
}
UInt_t AliITSOnlineSPDphysAnalyzer::GetNrDeadChips() {
  // returns nr of dead chips
  if (!fbDeadProcessed) ProcessDeadPixels();
  return fNrDeadChips;
}
UInt_t AliITSOnlineSPDphysAnalyzer::GetNrInefficientChips() {
  // returns nr of inefficient chips
  if (!fbDeadProcessed) ProcessDeadPixels();
  return fNrInefficientChips;
}
UInt_t AliITSOnlineSPDphysAnalyzer::GetNrNeedsMoreStatChips() {
  // returns nr of needs more stat chips
  if (!fbDeadProcessed) ProcessDeadPixels();
  return 60-fNrEnoughStatChips-fNrDeadChips-fNrInefficientChips;
}

UInt_t AliITSOnlineSPDphysAnalyzer::GetEqNr() const {
  // returns the eq nr of phys obj
  if (fPhysObj!=NULL) return fPhysObj->GetEqNr(); 
  else return 999;
}

UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEvents() const {
  // returns the nr of events of phys obj
  if (fPhysObj!=NULL) return fPhysObj->GetNrEvents(); 
  else return 0;
}

void AliITSOnlineSPDphysAnalyzer::Exponent(Double_t &val, Int_t &valExp) const {
  // put double in format with val and exp so that 1<val<10 - The actual value is val*10e(valExp)
  while (val>10) {
    val/=10;
    valExp++;
  }
  while (val<1) {
    val*=10;
    valExp--;
  }
}


TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapTot() {
  // creates and returns a pointer to a hitmap histo (half sector style a la spdmood)
  if (fPhysObj==NULL) {
    Error("AliITSOnlineSPDphysAnalyzer::GetHitMapTot","No data!");
    return NULL;
  }
  TString histoname = Form("Eq %d",GetEqNr());
  TH2F* fHitMapTot = new TH2F(histoname.Data(),histoname.Data(),32*10,-0.5,32*10-0.5,256*6,-0.5,256*6-0.5);
  fHitMapTot->SetNdivisions(-10,"X");
  fHitMapTot->SetNdivisions(-006,"Y");
  fHitMapTot->SetTickLength(0,"X");
  fHitMapTot->SetTickLength(0,"Y");
  fHitMapTot->GetXaxis()->SetLabelColor(gStyle->GetCanvasColor());
  fHitMapTot->GetYaxis()->SetLabelColor(gStyle->GetCanvasColor());
  for (UInt_t hs=0; hs<6; hs++) {
    for (UInt_t chipNr=0; chipNr<10; chipNr++) {
      for (UInt_t col=0; col<32; col++) {
	for (UInt_t row=0; row<256; row++) {
	  fHitMapTot->Fill(chipNr*32+col,(5-hs)*256+row,fPhysObj->GetHits(hs,chipNr,col,row));
	}
      }
    }
  }
  return fHitMapTot;
}

TH2F* AliITSOnlineSPDphysAnalyzer::GetHitMapChip(UInt_t hs, UInt_t chip) {
  // creates and returns a pointer to a hitmap histo (chip style a la spdmood)
  if (fPhysObj==NULL) {
    Error("AliITSOnlineSPDphysAnalyzer::GetHitMapChip","No data!");
    return NULL;
  }

  TString histoName;
  TString histoTitle;
  histoName = Form("fChipHisto_%d_%d_%d", GetEqNr(), hs, chip);
  histoTitle = Form("Eq ID %d, Half Stave %d, Chip %d", GetEqNr(), hs, chip);

  TH2F *returnHisto = new TH2F(histoName.Data(), histoTitle.Data(), 32, -0.5, 31.5, 256, -0.5, 255.5);
  returnHisto->SetMinimum(0);
  for (UInt_t col=0; col<32; col++) {
    for (UInt_t row=0; row<256; row++) {
      returnHisto->Fill(col,row,fPhysObj->GetHits(hs,chip,col,row));
    }
  }

  return returnHisto;
}
Commit	Line	Data
6727e2db	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(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
	42	fMinEventsForNoisy(10000),fMinEventsForDead(10000),
	43	fDefinitelyNoisyRatio(0.1),
	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(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
	55	fMinEventsForNoisy(10000),fMinEventsForDead(10000),
	56	fDefinitelyNoisyRatio(0.1),
	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(1),fThreshNoisyExp(-12),fThreshDead(1),fThreshDeadExp(-12),
68	fMinEventsForNoisy(10000),fMinEventsForDead(10000),
69	fDefinitelyNoisyRatio(0.1),
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	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;
83
84	Init();
85	}
86
87	AliITSOnlineSPDphysAnalyzer::~AliITSOnlineSPDphysAnalyzer() {
88	// destructor
89	if (fPhysObj!=NULL) delete fPhysObj;
90	}
91
92	AliITSOnlineSPDphysAnalyzer& AliITSOnlineSPDphysAnalyzer::operator=(const AliITSOnlineSPDphysAnalyzer& handle) {
93	// assignment operator, only copies the filename and params (not the processed data)
94	if (this!=&handle) {
95	if (fPhysObj!=NULL) delete fPhysObj;
96
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;
107
108	fPhysObj = NULL;
109	fHandler = NULL;
110	fNrEnoughStatChips = 0;
111	fNrDeadChips = 0;
112	fNrInefficientChips = 0;
113	fNrEqHits = 0;
114	fbDeadProcessed = kFALSE;
115
116	Init();
117	}
118	return *this;
119	}
120
121	void AliITSOnlineSPDphysAnalyzer::Init() {
122	// initialize container obj
123	FILE* fp0 = fopen(fFileName.Data(), "r");
124	if (fp0 == NULL) {
125	return;
126	}
127	else {
128	fclose(fp0);
129	}
130	fPhysObj = new AliITSOnlineSPDphys(fFileName.Data());
131	}
132
133	void AliITSOnlineSPDphysAnalyzer::SetParam(const Char_t pname, const Char_t pval) {
134	// set a parameter
135	TString name = pname;
136	TString val = pval;
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;
141	fThreshNoisyExp = 0;
142	Exponent(fThreshNoisy,fThreshNoisyExp);
143	fThreshNoisy/=(2061032256);
144	Exponent(fThreshNoisy,fThreshNoisyExp);
145	}
146	else if (name.CompareTo("MistakeProbabilityDead")==0) {
147	Double_t mistakeProbabilityDead = val.Atof();
148	fThreshDead = mistakeProbabilityDead;
149	fThreshDeadExp = 0;
150	Exponent(fThreshDead,fThreshDeadExp);
151	fThreshDead/=(2061032256);
152	Exponent(fThreshDead,fThreshDeadExp);
153	}
154	else if (name.CompareTo("fMinEventsForNoisy")==0) {
155	fMinEventsForNoisy = val.Atoi();
156	}
157	else if (name.CompareTo("fMinEventsForDead")==0) {
158	fMinEventsForDead = val.Atoi();
159	}
160	else if (name.CompareTo("fDefinitelyNoisyRatio")==0) {
161	fDefinitelyNoisyRatio = val.Atof();
162	}
163	else if (name.CompareTo("fMinNrEqHitsForDeadChips")==0) {
164	fMinNrEqHitsForDeadChips = val.Atof();
165	}
166	else if (name.CompareTo("fRatioToMeanForInefficientChip")==0) {
167	fRatioToMeanForInefficientChip = val.Atof();
168	}
169	else {
170	Error("AliITSOnlineSPDphysAnalyzer::SetParam","Parameter %s in configuration file unknown.",name.Data());
171	}
172	}
173
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);
177	ifstream paramsFile;
178	paramsFile.open(paramsFileName, ifstream::in);
179	if (paramsFile.fail()) {
180	printf("No config file (%s) present. Using default tuning parameters.\n",paramsFileName.Data());
181	}
182	else {
183	while(1) {
184	Char_t paramN[50];
185	Char_t paramV[50];
186	paramsFile >> paramN;
187	if (paramsFile.eof()) break;
188	paramsFile >> paramV;
189	SetParam(paramN,paramV);
190	if (paramsFile.eof()) break;
191	}
192	paramsFile.close();
193	}
194	}
195
196
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!");
201	return 0;
202	}
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);
206	return 0;
207	}
208	// handler should be initialized
209	if (fHandler==NULL) {
210	Error("AliITSOnlineSPDphysAnalyzer::ProcessNoisyPixels","Calibration handler is not initialized!");
211	return 0;
212	}
213
214	UInt_t nrEnoughStatChips = 0;
215
216	for (UInt_t hs=0; hs<6; hs++) {
217	for (UInt_t chip=0; chip<10; chip++) {
218
219	UInt_t nrPixels = 0;
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
227	nrPixels++;
228	if (nrHits>fDefinitelyNoisyRatio*GetNrEvents()) {
229	fHandler->SetNoisyPixel(GetEqNr(),hs,chip,col,row);
230	nrPixels--;
231	nrChipHits-=nrHits;
232	}
233	else {
234	if (nrMostHits<nrHits) nrMostHits=nrHits;
235	}
236	}
237	}
238
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) {
243
244	Double_t p = (Double_t)nrChipHits/nrPixels/n;
245
246	Double_t bin = 1;
247	Int_t binExp = 0;
248	// Bin(n,k=0):
249	for (UInt_t i=0; i<n; i++) {
250	bin*=(1-p);
251	Exponent(bin,binExp);
252	}
253	// Bin(n,k)
254	UInt_t k=1;
255	while (((binExp>fThreshNoisyExp \|\| (binExp==fThreshNoisyExp && bin>fThreshNoisy)) \|\| k<n*p) && k<=n) {
256	k++;
257	bin = bin(n-k+1)/kp/(1-p);
258	Exponent(bin,binExp);
259	}
260
261	// can we find noisy pixels...?
262	if (k<=n) {
263	// printf("eq %d , hs %d , chip %d : Noisy level = %d\n",GetEqNr(),hs,chip,k);
264	nrEnoughStatChips++;
265	// add noisy pixels to handler
266	UInt_t noiseLimit=k;
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);
273	}
274	}
275	}
276	}
277	}
278	}
279
280	}
281
282	} // for chip
283	} // for hs
284
285	return nrEnoughStatChips;
286	}
287
288
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!");
293	return 0;
294	}
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);
298	return 0;
299	}
300	// handler should be initialized
301	if (fHandler==NULL) {
302	Error("AliITSOnlineSPDphysAnalyzer::ProcessDeadPixels","Calibration handler is not initialized!");
303	return 0;
304	}
305
306	AliITSOnlineCalibrationSPDhandler *deadChipHandler = new AliITSOnlineCalibrationSPDhandler();
307	AliITSIntMap* possiblyIneff = new AliITSIntMap();
308
309	fNrEnoughStatChips = 0;
310	fNrDeadChips = 0;
311	fNrInefficientChips = 0;
312	UInt_t nrPossiblyDeadChips = 0;
313	fNrEqHits = 0;
314
315	for (UInt_t hs=0; hs<6; hs++) {
316	for (UInt_t chip=0; chip<10; chip++) {
317	Bool_t good=kFALSE;
318
319	AliITSOnlineCalibrationSPDhandler *deadPixelHandler = new AliITSOnlineCalibrationSPDhandler();
320
321	UInt_t nrPossiblyDeadPixels = 0;
322	UInt_t nrPixels = 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
330	nrPixels++;
331	if (nrHits==0) {
332	nrPossiblyDeadPixels++;
333	deadPixelHandler->SetDeadPixel(GetEqNr(),hs,chip,col,row);
334	}
335	}
336	else {
337	nrChipHits -= nrHits; // this is needed when running offline (online nrHits should be 0 already)
338	}
339	}
340	}
341	fNrEqHits+=nrChipHits;
342
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) {
347	if (nrChipHits>0) {
348	fHandler->UnSetDeadChip(GetEqNr(),hs,chip);
349	}
350	}
351	else {
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);
359	}
360	}
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));
364	}
365	delete tmpHand;
366	}
367
368
369
370
371	if (nrPossiblyDeadPixels==0) {
372	// no need to see if we have enough statistics...
373	fNrEnoughStatChips++;
374	good=kTRUE;
375	delete deadPixelHandler;
376	good=kTRUE;
377	printf("%3d",good);
378	if (chip==9) printf("\n");
379	continue;
380	}
381
382	if (nrChipHits==0) {
383	nrPossiblyDeadChips++;
384	deadChipHandler->SetDeadChip(GetEqNr(),hs,chip);
385	delete deadPixelHandler;
386	good=kFALSE;
387	printf("%3d",good);
388	if (chip==9) printf("\n");
389	continue;
390	}
391
392	// Binomial with n events and probability p for pixel hit
393	UInt_t n = GetNrEvents();
394	if (nrPixels>0 && n>0) {
395
396	Double_t p = (Double_t)nrChipHits/nrPixels/n;
397
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);
404	}
405
406	// can we find dead pixels...?
407	if (falselyDeadProbExp<fThreshDeadExp \|\| (falselyDeadProbExp==fThreshDeadExp && falselyDeadProb<fThreshDead)) {
408	fNrEnoughStatChips++;
409	good=kTRUE;
410	// add dead pixels to handler
411	fHandler->AddDeadFrom(deadPixelHandler);
412	}
413	else {
414	// this might be an inefficient chip
415	possiblyIneff->Insert(hs*10+chip,nrChipHits);
416	}
417
418	}
419	else {
420	if (n>0) {
421	// this is a completely noisy chip... put in category enough stat
422	fNrEnoughStatChips++;
423	good=kTRUE;
424	}
425	}
426
427	delete deadPixelHandler;
428
429	printf("%3d",good);
430	if (chip==9) printf("\n");
431
432	} // for chip
433	} // for hs
434
435
436	// dead chips?
437	if (fNrEqHits>fMinNrEqHitsForDeadChips) {
438	fHandler->AddDeadFrom(deadChipHandler);
439	fNrDeadChips=nrPossiblyDeadChips;
440	}
441	else {
442	fNrDeadChips=0;
443	}
444	delete deadChipHandler;
445
446
447	// inefficient chips?
448	Int_t key,val;
449	while (possiblyIneff->Pop(key,val)) {
450	if (val<fNrEqHits/60*fRatioToMeanForInefficientChip) {
451	fNrInefficientChips++;
452	}
453	}
454	delete possiblyIneff;
455
456
457	fbDeadProcessed = kTRUE;
458
459	return fNrEnoughStatChips;
460	}
461
462
463	UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEnoughStatChips() {
464	// returns nr of enough stat chips
465	if (!fbDeadProcessed) ProcessDeadPixels();
466	return fNrEnoughStatChips;
467	}
468	UInt_t AliITSOnlineSPDphysAnalyzer::GetNrDeadChips() {
469	// returns nr of dead chips
470	if (!fbDeadProcessed) ProcessDeadPixels();
471	return fNrDeadChips;
472	}
473	UInt_t AliITSOnlineSPDphysAnalyzer::GetNrInefficientChips() {
474	// returns nr of inefficient chips
475	if (!fbDeadProcessed) ProcessDeadPixels();
476	return fNrInefficientChips;
477	}
478	UInt_t AliITSOnlineSPDphysAnalyzer::GetNrNeedsMoreStatChips() {
479	// returns nr of needs more stat chips
480	if (!fbDeadProcessed) ProcessDeadPixels();
481	return 60-fNrEnoughStatChips-fNrDeadChips-fNrInefficientChips;
482	}
483
484	UInt_t AliITSOnlineSPDphysAnalyzer::GetEqNr() const {
485	// returns the eq nr of phys obj
486	if (fPhysObj!=NULL) return fPhysObj->GetEqNr();
487	else return 999;
488	}
489
490	UInt_t AliITSOnlineSPDphysAnalyzer::GetNrEvents() const {
491	// returns the nr of events of phys obj
492	if (fPhysObj!=NULL) return fPhysObj->GetNrEvents();
493	else return 0;
494	}
495
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)
498	while (val>10) {
499	val/=10;
500	valExp++;
501	}
502	while (val<1) {
503	val*=10;
504	valExp--;
505	}
506	}
507
508
509
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!");
514	return NULL;
515	}
516	TString histoname = Form("Eq %d",GetEqNr());
517	TH2F* fHitMapTot = new TH2F(histoname.Data(),histoname.Data(),3210,-0.5,3210-0.5,2566,-0.5,2566-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(chipNr32+col,(5-hs)256+row,fPhysObj->GetHits(hs,chipNr,col,row));
529	}
530	}
531	}
532	}
533	return fHitMapTot;
534	}
535
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!");
540	return NULL;
541	}
542
543	TString histoName;
544	TString histoTitle;
545	histoName = Form("fChipHisto_%d_%d_%d", GetEqNr(), hs, chip);
546	histoTitle = Form("Eq ID %d, Half Stave %d, Chip %d", GetEqNr(), hs, chip);
547
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));
553	}
554	}
555
556	return returnHisto;
557	}