end-of-line normalization
[u/mrichter/AliRoot.git] / TPC / Calib / AliTPCCalibKr.cxx
index 8ecdbbd..3b07e13 100644 (file)
-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-\r
-\r
-//Root includes\r
-#include <TH1F.h>\r
-#include <TH1D.h>\r
-#include <TH2F.h>\r
-#include <TH3F.h>\r
-#include <TString.h>\r
-#include <TMath.h>\r
-#include <TF1.h>\r
-#include <TRandom.h>\r
-#include <TDirectory.h>\r
-#include <TFile.h>\r
-#include <TAxis.h>\r
-//AliRoot includes\r
-#include "AliRawReader.h"\r
-#include "AliRawReaderRoot.h"\r
-#include "AliRawReaderDate.h"\r
-#include "AliTPCCalROC.h"\r
-#include "AliTPCCalPad.h"\r
-#include "AliTPCROC.h"\r
-#include "AliMathBase.h"\r
-#include "TTreeStream.h"\r
-\r
-//date\r
-#include "event.h"\r
-\r
-//header file\r
-#include "AliTPCCalibKr.h"\r
-\r
-//----------------------------------------------------------------------------\r
-// The AliTPCCalibKr class description (TPC Kr calibration).\r
-//\r
-//\r
-// The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),\r
-// its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).   \r
-// \r
-// The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.\r
-// These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration \r
-// parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.\r
-// In addition the debugCalibKr.root file with debug information is created.\r
-//\r
-\r
-/*\r
\r
-// Usage example:\r
-//\r
-\r
-// 1. Analyse output histograms:\r
-TFile f("outHistFile.root");\r
-AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr");\r
-obj->SetRadius(0,0);\r
-obj->SetStep(1,1);\r
-obj->Analyse();\r
-\r
-// 2. See calibration parameters e.g.:\r
-TFile f("calibKr.root");\r
-spectrMean->GetCalROC(70)->GetValue(40,40);\r
-fitMean->GetCalROC(70)->GetValue(40,40);\r
-\r
-// 3. See debug information e.g.:\r
-TFile f("debugCalibKr.root");\r
-.ls;\r
-\r
-// -- Print calibKr TTree content \r
-calibKr->Print();\r
-\r
-// -- Draw calibKr TTree variables\r
-calibKr.Draw("fitMean");\r
-\r
-*/\r
-\r
-\r
-//\r
-// Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)\r
-//-----------------------------------------------------------------------------\r
-\r
-ClassImp(AliTPCCalibKr)\r
-\r
-AliTPCCalibKr::AliTPCCalibKr() : \r
-  TObject(),\r
-  fASide(kTRUE),\r
-  fCSide(kTRUE),\r
-  fHistoKrArray(72),\r
-  fADCOverClustSizeMin(0.0),\r
-  fADCOverClustSizeMax(1.0e9),\r
-  fMaxADCOverClustADCMin(0.0),\r
-  fMaxADCOverClustADCMax(1.0e9),\r
-  fTimeMin(0.0),\r
-  fTimeMax(1.0e9),\r
-  fClustSizeMin(0.0),\r
-  fClustSizeMax(1.0e9),\r
-  fTimebinRmsIrocMin(0.0),\r
-  fPadRmsIrocMin(0.0),\r
-  fRowRmsIrocMin(0.0),\r
-  fClusterPadSize1DIrocMax(200),\r
-  fCurveCoefficientIroc(1.0e9),\r
-  fTimebinRmsOrocMin(0.0),\r
-  fPadRmsOrocMin(0.0),\r
-  fRowRmsOrocMin(0.0),\r
-  fClusterPadSize1DOrocMax(200),\r
-  fCurveCoefficientOroc(1.0e9),\r
-  fIrocHistogramMin(100.),\r
-  fIrocHistogramMax(6000.),\r
-  fIrocHistogramNbins(200),\r
-  fOrocHistogramMin(100.),\r
-  fOrocHistogramMax(5500.),\r
-  fOrocHistogramNbins(200),\r
-  fRowRadius(0),\r
-  fPadRadius(0),\r
-  fRowStep(1),\r
-  fPadStep(1)\r
-\r
-{\r
-  //\r
-  // default constructor\r
-  //\r
-\r
-  // TObjArray with histograms\r
-  fHistoKrArray.SetOwner(kTRUE); // is owner of histograms\r
-  fHistoKrArray.Clear();\r
\r
-  // init cuts (by Stefan)\r
-//  SetADCOverClustSizeRange(7,200);\r
-//  SetMaxADCOverClustADCRange(0.01,0.4);\r
-//  SetTimeRange(200,600);\r
-//  SetClustSizeRange(6,200);\r
-\r
-  //init  cuts (by Adam)\r
-  SetTimebinRmsMin(1.6,0.8);\r
-  SetPadRmsMin(0.825,0.55);\r
-  SetRowRmsMin(0.1,0.1);\r
-  SetClusterPadSize1DMax(15,11);\r
-  SetCurveCoefficient(1.,2.);\r
\r
-  //set histograms settings\r
-  SetIrocHistogram(200,100,6000);\r
-  SetOrocHistogram(200,100,5500);\r
-\r
-  // init histograms\r
-  //Init();\r
-}\r
-\r
-//_____________________________________________________________________\r
-AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) : \r
-  TObject(pad),\r
-  \r
-  fASide(pad.fASide),\r
-  fCSide(pad.fCSide),\r
-  fHistoKrArray(72),\r
-  fADCOverClustSizeMin(pad.fADCOverClustSizeMin),\r
-  fADCOverClustSizeMax(pad.fADCOverClustSizeMax),\r
-  fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),\r
-  fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),\r
-  fTimeMin(pad.fTimeMin),\r
-  fTimeMax(pad.fTimeMax),\r
-  fClustSizeMin(pad.fClustSizeMin),\r
-  fClustSizeMax(pad.fClustSizeMax),\r
-  fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),\r
-  fPadRmsIrocMin(pad.fPadRmsIrocMin),\r
-  fRowRmsIrocMin(pad.fRowRmsIrocMin),\r
-  fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),\r
-  fCurveCoefficientIroc(pad.fCurveCoefficientIroc),\r
-  fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),\r
-  fPadRmsOrocMin(pad.fPadRmsOrocMin),\r
-  fRowRmsOrocMin(pad.fRowRmsOrocMin),\r
-  fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),\r
-  fCurveCoefficientOroc(pad.fCurveCoefficientOroc),\r
-  fIrocHistogramMin(pad.fIrocHistogramMin),\r
-  fIrocHistogramMax(pad.fIrocHistogramMax),\r
-  fIrocHistogramNbins(pad.fIrocHistogramNbins),\r
-  fOrocHistogramMin(pad.fOrocHistogramMin),\r
-  fOrocHistogramMax(pad.fOrocHistogramMax),\r
-  fOrocHistogramNbins(pad.fOrocHistogramNbins),\r
-  fRowRadius(pad.fRowRadius),\r
-  fPadRadius(pad.fPadRadius),\r
-  fRowStep(pad.fRowStep),\r
-  fPadStep(pad.fPadStep)\r
-\r
-{\r
-  // copy constructor\r
\r
-  // TObjArray with histograms\r
-  fHistoKrArray.SetOwner(kTRUE); // is owner of histograms\r
-  fHistoKrArray.Clear();\r
-\r
-  for (Int_t iSec = 0; iSec < 72; ++iSec) \r
-  {\r
-    TH3F *hOld = pad.GetHistoKr(iSec);\r
-       if(hOld) {\r
-      TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) ); \r
-      fHistoKrArray.AddAt(hNew,iSec);\r
-       }\r
-  }\r
-}\r
-\r
-//_____________________________________________________________________\r
-AliTPCCalibKr::~AliTPCCalibKr() \r
-{\r
-  //\r
-  // destructor\r
-  //\r
-\r
-  // for (Int_t iSec = 0; iSec < 72; ++iSec) {\r
-  //     if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);\r
-  // }\r
-  fHistoKrArray.Delete();\r
-}\r
-\r
-//_____________________________________________________________________\r
-AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)\r
-{\r
-  // assignment operator\r
-\r
-  if (&source == this) return *this;\r
-  new (this) AliTPCCalibKr(source);\r
-\r
-  return *this;\r
-}\r
-\r
-//_____________________________________________________________________\r
-void AliTPCCalibKr::Init()\r
-{\r
-  // \r
-  // init output histograms \r
-  //\r
-\r
-  // add histograms to the TObjArray\r
-  for(Int_t i=0; i<72; ++i) {\r
-    \r
-    // C - side\r
-    if(IsCSide(i) == kTRUE && fCSide == kTRUE) {\r
-      TH3F *hist = CreateHisto(i);\r
-      if(hist) fHistoKrArray.AddAt(hist,i);\r
-    }\r
-    \r
-    // A - side\r
-    if(IsCSide(i) == kFALSE && fASide == kTRUE) {\r
-      TH3F *hist = CreateHisto(i);\r
-      if(hist) fHistoKrArray.AddAt(hist,i);\r
-    }\r
-  }\r
-}\r
\r
-//_____________________________________________________________________\r
-Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)\r
-{\r
-  //\r
-  // process events \r
-  // call event by event\r
-  //\r
-\r
-  if(cluster) Update(cluster);\r
-  else return kFALSE;\r
\r
- return kTRUE;\r
-}\r
-\r
-//_____________________________________________________________________\r
-TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)\r
-{\r
-    //\r
-    // create new histogram\r
-       //\r
-    char name[256];\r
-       TH3F *h;\r
-\r
-       snprintf(name,256,"ADCcluster_ch%d",chamber);\r
-\r
-    if( IsIROC(chamber) == kTRUE ) \r
-       {\r
-         h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);\r
-       } else {\r
-         h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);\r
-       }\r
-    h->SetXTitle("padrow");\r
-    h->SetYTitle("pad");\r
-    h->SetZTitle("fADC");\r
-\r
-return h;\r
-}\r
-\r
-//_____________________________________________________________________\r
-Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)\r
-{\r
-// check if IROCs\r
-// returns kTRUE if IROCs and kFALSE if OROCs \r
-\r
-  if(chamber>=0 && chamber<36) return kTRUE;\r
-\r
-return kFALSE;\r
-}\r
-\r
-//_____________________________________________________________________\r
-Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)\r
-{\r
-// check if C side\r
-// returns kTRUE if C side and kFALSE if A side\r
-\r
-  if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;\r
-\r
-return kFALSE;\r
-}\r
\r
-//_____________________________________________________________________\r
-Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)\r
-{\r
-  //\r
-  // fill existing histograms\r
-  //\r
-\r
-  if (!Accept(cl)) return kFALSE;\r
-  TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());\r
-  if(!h) return kFALSE;\r
-  \r
-  h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());\r
-  \r
-  return kTRUE;\r
-}\r
-\r
-//_____________________________________________________________________\r
-Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){\r
-  //\r
-  // cuts\r
-  //\r
-  /*\r
-    TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings - \r
-    TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal\r
-    TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal\r
-    TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal\r
-    TCut cutR4("cutR4","fMax.fTime>200");   // noise removal\r
-    TCut cutR5("cutR5","fMax.fTime<600");   // noise removal\r
-    TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks\r
-    TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;\r
-  */\r
-  /*\r
-  //R0\r
-  if ((float)cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;\r
-  // R1\r
-  if ((float)cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;\r
-  //R2\r
-  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;\r
-  //R3\r
-  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;\r
-  //R4\r
-  if (cl->GetMax().GetTime() < 200) return kFALSE;\r
-  //R5\r
-  if (cl->GetMax().GetTime() > 600) return kFALSE;\r
-  //S1\r
-  if (cl->GetSize()>200) return kFALSE;\r
-  if (cl->GetSize()<6)  return kFALSE;\r
-\r
-  SetADCOverClustSizeRange(7,200);\r
-  SetMaxADCOverClustADCRange(0.01,0.4);\r
-  SetTimeRange(200,600);\r
-  SetClustSizeRange(6,200);\r
-  */\r
-\r
-  //R0\r
-  if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax)        return kFALSE;\r
-  // R1\r
-  if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin)          return kFALSE;\r
-  //R2\r
-  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax)  return kFALSE;\r
-  //R3\r
-  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;\r
-  //R4\r
-  if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;\r
-  //R5\r
-  if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;\r
-  //S1\r
-  if (cl->GetSize()>fClustSizeMax) return kFALSE;\r
-  if (cl->GetSize()<fClustSizeMin)  return kFALSE;\r
-\r
-  //\r
-  // cuts by Adam\r
-  //\r
-  /*\r
-    TCut cutAI0("cutAI0","fTimebinRMS>1.6");\r
-    TCut cutAI1("cutAI1","fPadRMS>0.825");  \r
-    TCut cutAI2("cutAI2","fRowRMS>0.1");    \r
-    TCut cutAI3("cutAI3","fPads1D<15");  \r
-    TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0"); \r
-\r
-    TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;\r
-\r
-    TCut cutAO0("cutAO0","fTimebinRMS>0.8");\r
-    TCut cutAO1("cutAO1","fPadRMS>0.55");  \r
-    TCut cutAO2("cutAO2","fRowRMS>0.1");    \r
-    TCut cutAO3("cutAO3","fPads1D<11");  \r
-    TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0"); \r
-\r
-    TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;\r
-    TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");\r
-\r
-  */\r
-  if(cl->GetSec()<36){  //IROCs\r
-    //AI0\r
-    if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;\r
-    //AI1\r
-    if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;\r
-    //AI2\r
-    if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;\r
-    //AI3\r
-    if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;\r
-    //AI4\r
-    if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;\r
-  }else{  //OROCs\r
-    //AO0\r
-    if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;\r
-    //AO1\r
-    if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;\r
-    //AO2\r
-    if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;\r
-    //AO3\r
-    if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;\r
-    //AO4\r
-    if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;\r
-  }\r
-\r
-  return kTRUE;\r
-\r
-}\r
-\r
-//_____________________________________________________________________\r
-TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const\r
-{\r
-  // get histograms from fHistoKrArray\r
-  return (TH3F*) fHistoKrArray.At(chamber);\r
-}\r
\r
-//_____________________________________________________________________\r
-void AliTPCCalibKr::Analyse() \r
-{\r
-  //\r
-  // analyse the histograms and extract krypton calibration parameters\r
-  //\r
-\r
-  // AliTPCCalPads that will contain the calibration parameters\r
-  AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");\r
-  AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");\r
-  AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");\r
-  AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");\r
-  AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");\r
-  AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");\r
-\r
-  // file stream for debugging purposes\r
-  TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");\r
-\r
-  // if entries in spectrum less than minEntries, then the fit won't be performed\r
-  Int_t minEntries = 1; //300;\r
-\r
-  Double_t windowFrac = 0.12;\r
-  // the 3d histogram will be projected on the pads given by the following window size\r
-  // set the numbers to 0 if you want to do a pad-by-pad calibration\r
-  UInt_t rowRadius = fRowRadius;//4;\r
-  UInt_t padRadius = fPadRadius;//4;\r
-  \r
-  // the step size by which pad and row are incremented is given by the following numbers\r
-  // set them to 1 if you want the finest granularity\r
-  UInt_t rowStep = fRowStep;//1;//2    // formerly: 2*rowRadius\r
-  UInt_t padStep = fPadStep;//1;//4    // formerly: 2*padRadius\r
-\r
-  for (Int_t chamber = 0; chamber < 72; chamber++) {\r
-    //if (chamber != 71) continue;\r
-    AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()\r
-    \r
-    // Usually I would traverse each pad, take the spectrum for its neighbourhood and\r
-    // obtain the calibration parameters. This takes very long, so instead I assign the same\r
-    // calibration values to the whole neighbourhood and then go on to the next neighbourhood.\r
-    UInt_t nRows = roc.GetNrows();\r
-    for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {\r
-      UInt_t nPads = roc.GetNPads(iRow);\r
-      //if (iRow >= 10) break;\r
-      for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {\r
-        //if (iPad >= 20) break;\r
-        TH3F* h = GetHistoKr(chamber);\r
-        if (!h) continue;\r
-        \r
-        // the 3d histogram will be projected on the pads given by the following bounds\r
-        // for rows and pads\r
-        Int_t rowLow = iRow - rowRadius;\r
-        UInt_t rowUp = iRow + rowRadius + rowStep-1;\r
-        Int_t padLow = iPad - padRadius;\r
-        UInt_t padUp = iPad + padRadius + padStep-1;\r
-        // if window goes out of chamber\r
-        if (rowLow < 0) rowLow = 0;\r
-        if (rowUp >= nRows) rowUp = nRows - 1;\r
-        if (padLow < 0) padLow = 0;\r
-        if (padUp >= nPads) padUp = nPads - 1;\r
-\r
-        // project the histogram\r
-        //TH1D* projH = h->ProjectionZ("projH", rowLow+1, rowUp+1, padLow+1, padUp+1); // SLOW\r
-        TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);\r
-    \r
-        // get the number of entries in the spectrum\r
-        Double_t entries = projH->GetEntries();\r
-        if (entries < minEntries) { delete projH; continue; }\r
-        \r
-        // get the two calibration parameters mean of spectrum and RMS of spectrum\r
-        Double_t histMean = projH->GetMean();\r
-        Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;\r
-    \r
-        // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted\r
-        Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());\r
-               Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);\r
-               Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);\r
-        Double_t integCharge =  projH->Integral(minBin,maxBin); \r
-\r
-        Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);\r
-\r
-        if (fitResult != 0) {\r
-          Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);\r
-          //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f,  %f\n", chamber, iRow, iPad, entries, maxEntries);\r
-    \r
-          delete projH;\r
-          continue;\r
-        }\r
-    \r
-        // get the two calibration parameters mean of gauss fit and sigma of gauss fit\r
-        TF1* gausFit = projH->GetFunction("gaus");\r
-        Double_t fitMean = gausFit->GetParameter(1);\r
-        Double_t fitRMS = gausFit->GetParameter(2);\r
-        Int_t numberFitPoints = gausFit->GetNumberFitPoints();\r
-        if (numberFitPoints == 0) continue;\r
-        Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;\r
-        delete gausFit;\r
-        delete projH;\r
-        if (fitMean <= 0) continue;\r
-        //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);\r
-    \r
-        // write the calibration parameters for each pad that the 3d histogram was projected onto\r
-        // (with considering the step size) to the CalPads\r
-        // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions\r
-        // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction\r
-        for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {\r
-          if (r < 0 || r >= (Int_t)nRows) continue;\r
-          UInt_t nPadsR = roc.GetNPads(r);\r
-          for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {\r
-            if (p < 0 || p >= (Int_t)nPadsR) continue;\r
-            spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);\r
-            spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);\r
-            fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);\r
-            fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);\r
-            fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);\r
-            entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);\r
-\r
-            (*debugStream) << "calibKr" <<\r
-              "sector=" << chamber <<          // chamber number\r
-              "row=" << r <<                   // row number\r
-              "pad=" << p <<                   // pad number\r
-              "histMean=" << histMean <<       // mean of the spectrum\r
-              "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean\r
-              "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak\r
-              "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak\r
-              "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit\r
-              "entries=" << entries <<         // number of entries for the spectrum\r
-              "\n";\r
-          }\r
-        }\r
-      }\r
-    }\r
-  }\r
-\r
-  TFile f("calibKr.root", "recreate");\r
-  spectrMeanCalPad->Write();\r
-  spectrRMSCalPad->Write();\r
-  fitMeanCalPad->Write();\r
-  fitRMSCalPad->Write();\r
-  fitNormChi2CalPad->Write();\r
-  entriesCalPad->Write();\r
-  f.Close();\r
-  delete spectrMeanCalPad;\r
-  delete spectrRMSCalPad;\r
-  delete fitMeanCalPad;\r
-  delete fitRMSCalPad;\r
-  delete fitNormChi2CalPad;\r
-  delete entriesCalPad;\r
-  delete debugStream;\r
-}\r
-\r
-//_____________________________________________________________________\r
-TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)\r
-{\r
-  // project the z-axis of a 3d histo to a specific range of the x- and y-axes,\r
-  // replaces TH3F::ProjectZ() to gain more speed\r
-\r
-  TAxis* xAxis = histo3D->GetXaxis();\r
-  TAxis* yAxis = histo3D->GetYaxis();\r
-  TAxis* zAxis = histo3D->GetZaxis();\r
-  Double_t zMinVal = zAxis->GetXmin();\r
-  Double_t zMaxVal = zAxis->GetXmax();\r
-  \r
-  Int_t nBinsZ = zAxis->GetNbins();\r
-  TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);\r
-\r
-  Int_t nx = xAxis->GetNbins()+2;\r
-  Int_t ny = yAxis->GetNbins()+2;\r
-  Int_t bin = 0;\r
-  Double_t entries = 0.;\r
-  for (Int_t x = xMin; x <= xMax; x++) {\r
-    for (Int_t y = yMin; y <= yMax; y++) {\r
-      for (Int_t z = 0; z <= nBinsZ+1; z++) {\r
-        bin = x + nx * (y + ny * z);\r
-        Double_t val = histo3D->GetBinContent(bin);\r
-        projH->Fill(zAxis->GetBinCenter(z), val);\r
-        entries += val;\r
-      }\r
-    }\r
-  }\r
-  projH->SetEntries((Long64_t)entries);\r
-  return projH;\r
-}\r
-\r
-//_____________________________________________________________________\r
-Long64_t AliTPCCalibKr::Merge(TCollection* list) {\r
-// merge function \r
-//\r
-\r
-if (!list)\r
-return 0;\r
-\r
-if (list->IsEmpty())\r
-return 1;\r
-\r
-TIterator* iter = list->MakeIterator();\r
-TObject* obj = 0;\r
-\r
-    iter->Reset();\r
-    Int_t count=0;\r
-       while((obj = iter->Next()) != 0)\r
-       {\r
-         AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);\r
-         if (entry == 0) continue; \r
-\r
-               for(int i=0; i<72; ++i) { \r
-                 if(IsCSide(i) == kTRUE && fCSide == kTRUE) { \r
-                 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  \r
-                 } \r
-\r
-                 if(IsCSide(i) == kFALSE && fASide == kTRUE) { \r
-                   ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  \r
-                 } \r
-               } \r
-\r
-         count++;\r
-       }\r
-\r
-return count;\r
-}\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, 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.                  *
+ **************************************************************************/
+
+
+//Root includes
+#include <TH1F.h>
+#include <TH1D.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TF1.h>
+#include <TRandom.h>
+#include <TDirectory.h>
+#include <TFile.h>
+#include <TAxis.h>
+//AliRoot includes
+#include "AliRawReader.h"
+#include "AliRawReaderRoot.h"
+#include "AliRawReaderDate.h"
+#include "AliTPCCalROC.h"
+#include "AliTPCCalPad.h"
+#include "AliTPCROC.h"
+#include "AliMathBase.h"
+#include "TTreeStream.h"
+
+//date
+#include "event.h"
+
+//header file
+#include "AliTPCCalibKr.h"
+
+//----------------------------------------------------------------------------
+// The AliTPCCalibKr class description (TPC Kr calibration).
+//
+//
+// The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
+// its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).   
+// 
+// The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
+// These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration 
+// parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
+// In addition the debugCalibKr.root file with debug information is created.
+//
+
+/*
+// Usage example:
+//
+
+// 1. Analyse output histograms:
+TFile f("outHistFile.root");
+AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr");
+obj->SetRadius(0,0);
+obj->SetStep(1,1);
+obj->Analyse();
+
+// 2. See calibration parameters e.g.:
+TFile f("calibKr.root");
+spectrMean->GetCalROC(70)->GetValue(40,40);
+fitMean->GetCalROC(70)->GetValue(40,40);
+
+// 3. See debug information e.g.:
+TFile f("debugCalibKr.root");
+.ls;
+
+// -- Print calibKr TTree content 
+calibKr->Print();
+
+// -- Draw calibKr TTree variables
+calibKr.Draw("fitMean");
+
+*/
+
+
+//
+// Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
+//-----------------------------------------------------------------------------
+
+ClassImp(AliTPCCalibKr)
+
+AliTPCCalibKr::AliTPCCalibKr() : 
+  TObject(),
+  fASide(kTRUE),
+  fCSide(kTRUE),
+  fHistoKrArray(72),
+  fADCOverClustSizeMin(0.0),
+  fADCOverClustSizeMax(1.0e9),
+  fMaxADCOverClustADCMin(0.0),
+  fMaxADCOverClustADCMax(1.0e9),
+  fTimeMin(0.0),
+  fTimeMax(1.0e9),
+  fClustSizeMin(0.0),
+  fClustSizeMax(1.0e9),
+  fTimebinRmsIrocMin(0.0),
+  fPadRmsIrocMin(0.0),
+  fRowRmsIrocMin(0.0),
+  fClusterPadSize1DIrocMax(200),
+  fCurveCoefficientIroc(1.0e9),
+  fTimebinRmsOrocMin(0.0),
+  fPadRmsOrocMin(0.0),
+  fRowRmsOrocMin(0.0),
+  fClusterPadSize1DOrocMax(200),
+  fCurveCoefficientOroc(1.0e9),
+  fIrocHistogramMin(100.),
+  fIrocHistogramMax(6000.),
+  fIrocHistogramNbins(200),
+  fOrocHistogramMin(100.),
+  fOrocHistogramMax(5500.),
+  fOrocHistogramNbins(200),
+  fRowRadius(0),
+  fPadRadius(0),
+  fRowStep(1),
+  fPadStep(1)
+
+{
+  //
+  // default constructor
+  //
+
+  // TObjArray with histograms
+  fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
+  fHistoKrArray.Clear();
+  // init cuts (by Stefan)
+//  SetADCOverClustSizeRange(7,200);
+//  SetMaxADCOverClustADCRange(0.01,0.4);
+//  SetTimeRange(200,600);
+//  SetClustSizeRange(6,200);
+
+  //init  cuts (by Adam)
+  SetTimebinRmsMin(1.6,0.8);
+  SetPadRmsMin(0.825,0.55);
+  SetRowRmsMin(0.1,0.1);
+  SetClusterPadSize1DMax(15,11);
+  SetCurveCoefficient(1.,2.);
+  //set histograms settings
+  SetIrocHistogram(200,100,6000);
+  SetOrocHistogram(200,100,5500);
+
+  // init histograms
+  //Init();
+}
+
+//_____________________________________________________________________
+AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) : 
+  TObject(pad),
+  
+  fASide(pad.fASide),
+  fCSide(pad.fCSide),
+  fHistoKrArray(72),
+  fADCOverClustSizeMin(pad.fADCOverClustSizeMin),
+  fADCOverClustSizeMax(pad.fADCOverClustSizeMax),
+  fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),
+  fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),
+  fTimeMin(pad.fTimeMin),
+  fTimeMax(pad.fTimeMax),
+  fClustSizeMin(pad.fClustSizeMin),
+  fClustSizeMax(pad.fClustSizeMax),
+  fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),
+  fPadRmsIrocMin(pad.fPadRmsIrocMin),
+  fRowRmsIrocMin(pad.fRowRmsIrocMin),
+  fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),
+  fCurveCoefficientIroc(pad.fCurveCoefficientIroc),
+  fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),
+  fPadRmsOrocMin(pad.fPadRmsOrocMin),
+  fRowRmsOrocMin(pad.fRowRmsOrocMin),
+  fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),
+  fCurveCoefficientOroc(pad.fCurveCoefficientOroc),
+  fIrocHistogramMin(pad.fIrocHistogramMin),
+  fIrocHistogramMax(pad.fIrocHistogramMax),
+  fIrocHistogramNbins(pad.fIrocHistogramNbins),
+  fOrocHistogramMin(pad.fOrocHistogramMin),
+  fOrocHistogramMax(pad.fOrocHistogramMax),
+  fOrocHistogramNbins(pad.fOrocHistogramNbins),
+  fRowRadius(pad.fRowRadius),
+  fPadRadius(pad.fPadRadius),
+  fRowStep(pad.fRowStep),
+  fPadStep(pad.fPadStep)
+
+{
+  // copy constructor
+  // TObjArray with histograms
+  fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
+  fHistoKrArray.Clear();
+
+  for (Int_t iSec = 0; iSec < 72; ++iSec) 
+  {
+    TH3F *hOld = pad.GetHistoKr(iSec);
+       if(hOld) {
+      TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) ); 
+      fHistoKrArray.AddAt(hNew,iSec);
+       }
+  }
+}
+
+//_____________________________________________________________________
+AliTPCCalibKr::~AliTPCCalibKr() 
+{
+  //
+  // destructor
+  //
+
+  // for (Int_t iSec = 0; iSec < 72; ++iSec) {
+  //     if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
+  // }
+  fHistoKrArray.Delete();
+}
+
+//_____________________________________________________________________
+AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)
+{
+  // assignment operator
+
+  if (&source == this) return *this;
+  new (this) AliTPCCalibKr(source);
+
+  return *this;
+}
+
+//_____________________________________________________________________
+void AliTPCCalibKr::Init()
+{
+  // 
+  // init output histograms 
+  //
+
+  // add histograms to the TObjArray
+  for(Int_t i=0; i<72; ++i) {
+    
+    // C - side
+    if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
+      TH3F *hist = CreateHisto(i);
+      if(hist) fHistoKrArray.AddAt(hist,i);
+    }
+    
+    // A - side
+    if(IsCSide(i) == kFALSE && fASide == kTRUE) {
+      TH3F *hist = CreateHisto(i);
+      if(hist) fHistoKrArray.AddAt(hist,i);
+    }
+  }
+}
+//_____________________________________________________________________
+Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
+{
+  //
+  // process events 
+  // call event by event
+  //
+
+  if(cluster) Update(cluster);
+  else return kFALSE;
+ return kTRUE;
+}
+
+//_____________________________________________________________________
+TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
+{
+    //
+    // create new histogram
+       //
+    char name[256];
+       TH3F *h;
+
+       snprintf(name,256,"ADCcluster_ch%d",chamber);
+
+    if( IsIROC(chamber) == kTRUE ) 
+       {
+         h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);
+       } else {
+         h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);
+       }
+    h->SetXTitle("padrow");
+    h->SetYTitle("pad");
+    h->SetZTitle("fADC");
+
+return h;
+}
+
+//_____________________________________________________________________
+Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
+{
+// check if IROCs
+// returns kTRUE if IROCs and kFALSE if OROCs 
+
+  if(chamber>=0 && chamber<36) return kTRUE;
+
+return kFALSE;
+}
+
+//_____________________________________________________________________
+Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
+{
+// check if C side
+// returns kTRUE if C side and kFALSE if A side
+
+  if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
+
+return kFALSE;
+}
+//_____________________________________________________________________
+Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)
+{
+  //
+  // fill existing histograms
+  //
+
+  if (!Accept(cl)) return kFALSE;
+  TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
+  if(!h) return kFALSE;
+  
+  h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
+  
+  return kTRUE;
+}
+
+//_____________________________________________________________________
+Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){
+  //
+  // cuts
+  //
+  /*
+    TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings - 
+    TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal
+    TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal
+    TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal
+    TCut cutR4("cutR4","fMax.fTime>200");   // noise removal
+    TCut cutR5("cutR5","fMax.fTime<600");   // noise removal
+    TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks
+    TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;
+  */
+  /*
+  //R0
+  if ((float)cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;
+  // R1
+  if ((float)cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;
+  //R2
+  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;
+  //R3
+  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
+  //R4
+  if (cl->GetMax().GetTime() < 200) return kFALSE;
+  //R5
+  if (cl->GetMax().GetTime() > 600) return kFALSE;
+  //S1
+  if (cl->GetSize()>200) return kFALSE;
+  if (cl->GetSize()<6)  return kFALSE;
+
+  SetADCOverClustSizeRange(7,200);
+  SetMaxADCOverClustADCRange(0.01,0.4);
+  SetTimeRange(200,600);
+  SetClustSizeRange(6,200);
+  */
+
+  //R0
+  if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax)        return kFALSE;
+  // R1
+  if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin)          return kFALSE;
+  //R2
+  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax)  return kFALSE;
+  //R3
+  if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
+  //R4
+  if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
+  //R5
+  if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
+  //S1
+  if (cl->GetSize()>fClustSizeMax) return kFALSE;
+  if (cl->GetSize()<fClustSizeMin)  return kFALSE;
+
+  //
+  // cuts by Adam
+  //
+  /*
+    TCut cutAI0("cutAI0","fTimebinRMS>1.6");
+    TCut cutAI1("cutAI1","fPadRMS>0.825");  
+    TCut cutAI2("cutAI2","fRowRMS>0.1");    
+    TCut cutAI3("cutAI3","fPads1D<15");  
+    TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0"); 
+
+    TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;
+
+    TCut cutAO0("cutAO0","fTimebinRMS>0.8");
+    TCut cutAO1("cutAO1","fPadRMS>0.55");  
+    TCut cutAO2("cutAO2","fRowRMS>0.1");    
+    TCut cutAO3("cutAO3","fPads1D<11");  
+    TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0"); 
+
+    TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;
+    TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");
+
+  */
+  if(cl->GetSec()<36){  //IROCs
+    //AI0
+    if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;
+    //AI1
+    if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;
+    //AI2
+    if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;
+    //AI3
+    if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;
+    //AI4
+    if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
+  }else{  //OROCs
+    //AO0
+    if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;
+    //AO1
+    if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;
+    //AO2
+    if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;
+    //AO3
+    if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;
+    //AO4
+    if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
+  }
+
+  return kTRUE;
+
+}
+
+//_____________________________________________________________________
+TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
+{
+  // get histograms from fHistoKrArray
+  return (TH3F*) fHistoKrArray.At(chamber);
+}
+//_____________________________________________________________________
+void AliTPCCalibKr::Analyse() 
+{
+  //
+  // analyse the histograms and extract krypton calibration parameters
+  //
+
+  // AliTPCCalPads that will contain the calibration parameters
+  AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
+  AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
+  AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
+  AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
+  AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
+  AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
+
+  // file stream for debugging purposes
+  TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
+
+  // if entries in spectrum less than minEntries, then the fit won't be performed
+  Int_t minEntries = 1; //300;
+
+  Double_t windowFrac = 0.12;
+  // the 3d histogram will be projected on the pads given by the following window size
+  // set the numbers to 0 if you want to do a pad-by-pad calibration
+  UInt_t rowRadius = fRowRadius;//4;
+  UInt_t padRadius = fPadRadius;//4;
+  
+  // the step size by which pad and row are incremented is given by the following numbers
+  // set them to 1 if you want the finest granularity
+  UInt_t rowStep = fRowStep;//1;//2    // formerly: 2*rowRadius
+  UInt_t padStep = fPadStep;//1;//4    // formerly: 2*padRadius
+
+  for (Int_t chamber = 0; chamber < 72; chamber++) {
+    //if (chamber != 71) continue;
+    AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()
+    
+    // Usually I would traverse each pad, take the spectrum for its neighbourhood and
+    // obtain the calibration parameters. This takes very long, so instead I assign the same
+    // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
+    UInt_t nRows = roc.GetNrows();
+    for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
+      UInt_t nPads = roc.GetNPads(iRow);
+      //if (iRow >= 10) break;
+      for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
+        //if (iPad >= 20) break;
+        TH3F* h = GetHistoKr(chamber);
+        if (!h) continue;
+        
+        // the 3d histogram will be projected on the pads given by the following bounds
+        // for rows and pads
+        Int_t rowLow = iRow - rowRadius;
+        UInt_t rowUp = iRow + rowRadius + rowStep-1;
+        Int_t padLow = iPad - padRadius;
+        UInt_t padUp = iPad + padRadius + padStep-1;
+        // if window goes out of chamber
+        if (rowLow < 0) rowLow = 0;
+        if (rowUp >= nRows) rowUp = nRows - 1;
+        if (padLow < 0) padLow = 0;
+        if (padUp >= nPads) padUp = nPads - 1;
+
+        // project the histogram
+        //TH1D* projH = h->ProjectionZ("projH", rowLow+1, rowUp+1, padLow+1, padUp+1); // SLOW
+        TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);
+    
+        // get the number of entries in the spectrum
+        Double_t entries = projH->GetEntries();
+        if (entries < minEntries) { delete projH; continue; }
+        
+        // get the two calibration parameters mean of spectrum and RMS of spectrum
+        Double_t histMean = projH->GetMean();
+        Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
+    
+        // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
+        Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
+               Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);
+               Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
+        Double_t integCharge =  projH->Integral(minBin,maxBin); 
+
+        Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
+
+        if (fitResult != 0) {
+          Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
+          //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f,  %f\n", chamber, iRow, iPad, entries, maxEntries);
+    
+          delete projH;
+          continue;
+        }
+    
+        // get the two calibration parameters mean of gauss fit and sigma of gauss fit
+        TF1* gausFit = projH->GetFunction("gaus");
+        Double_t fitMean = gausFit->GetParameter(1);
+        Double_t fitRMS = gausFit->GetParameter(2);
+        Int_t numberFitPoints = gausFit->GetNumberFitPoints();
+        if (numberFitPoints == 0) continue;
+        Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
+        delete gausFit;
+        delete projH;
+        if (fitMean <= 0) continue;
+        //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
+    
+        // write the calibration parameters for each pad that the 3d histogram was projected onto
+        // (with considering the step size) to the CalPads
+        // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
+        // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
+        for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
+          if (r < 0 || r >= (Int_t)nRows) continue;
+          UInt_t nPadsR = roc.GetNPads(r);
+          for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
+            if (p < 0 || p >= (Int_t)nPadsR) continue;
+            spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
+            spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
+            fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
+            fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
+            fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
+            entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
+
+            (*debugStream) << "calibKr" <<
+              "sector=" << chamber <<          // chamber number
+              "row=" << r <<                   // row number
+              "pad=" << p <<                   // pad number
+              "histMean=" << histMean <<       // mean of the spectrum
+              "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean
+              "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak
+              "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak
+              "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit
+              "entries=" << entries <<         // number of entries for the spectrum
+              "\n";
+          }
+        }
+      }
+    }
+  }
+
+  TFile f("calibKr.root", "recreate");
+  spectrMeanCalPad->Write();
+  spectrRMSCalPad->Write();
+  fitMeanCalPad->Write();
+  fitRMSCalPad->Write();
+  fitNormChi2CalPad->Write();
+  entriesCalPad->Write();
+  f.Close();
+  delete spectrMeanCalPad;
+  delete spectrRMSCalPad;
+  delete fitMeanCalPad;
+  delete fitRMSCalPad;
+  delete fitNormChi2CalPad;
+  delete entriesCalPad;
+  delete debugStream;
+}
+
+//_____________________________________________________________________
+TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
+{
+  // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
+  // replaces TH3F::ProjectZ() to gain more speed
+
+  TAxis* xAxis = histo3D->GetXaxis();
+  TAxis* yAxis = histo3D->GetYaxis();
+  TAxis* zAxis = histo3D->GetZaxis();
+  Double_t zMinVal = zAxis->GetXmin();
+  Double_t zMaxVal = zAxis->GetXmax();
+  
+  Int_t nBinsZ = zAxis->GetNbins();
+  TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
+
+  Int_t nx = xAxis->GetNbins()+2;
+  Int_t ny = yAxis->GetNbins()+2;
+  Int_t bin = 0;
+  Double_t entries = 0.;
+  for (Int_t x = xMin; x <= xMax; x++) {
+    for (Int_t y = yMin; y <= yMax; y++) {
+      for (Int_t z = 0; z <= nBinsZ+1; z++) {
+        bin = x + nx * (y + ny * z);
+        Double_t val = histo3D->GetBinContent(bin);
+        projH->Fill(zAxis->GetBinCenter(z), val);
+        entries += val;
+      }
+    }
+  }
+  projH->SetEntries((Long64_t)entries);
+  return projH;
+}
+
+//_____________________________________________________________________
+Long64_t AliTPCCalibKr::Merge(TCollection* list) {
+// merge function 
+//
+
+if (!list)
+return 0;
+
+if (list->IsEmpty())
+return 1;
+
+TIterator* iter = list->MakeIterator();
+TObject* obj = 0;
+
+    iter->Reset();
+    Int_t count=0;
+       while((obj = iter->Next()) != 0)
+       {
+         AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
+         if (entry == 0) continue; 
+
+               for(int i=0; i<72; ++i) { 
+                 if(IsCSide(i) == kTRUE && fCSide == kTRUE) { 
+                 ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  
+                 } 
+
+                 if(IsCSide(i) == kFALSE && fASide == kTRUE) { 
+                   ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );  
+                 } 
+               } 
+
+         count++;
+       }
+
+return count;
+}