TPC/TPCcalib/AliTPCCalibKr.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //Root includes
  18 #include <TH1F.h>
  19 #include <TH1D.h>
  20 #include <TH2F.h>
  21 #include <TH3F.h>
  22 #include <TString.h>
  23 #include <TMath.h>
  24 #include <TF1.h>
  25 #include <TRandom.h>
  26 #include <TDirectory.h>
  27 #include <TFile.h>
  28 #include <TAxis.h>
  29 //AliRoot includes
  30 #include "AliRawReader.h"
  31 #include "AliRawReaderRoot.h"
  32 #include "AliRawReaderDate.h"
  33 #include "AliTPCCalROC.h"
  34 #include "AliTPCCalPad.h"
  35 #include "AliTPCROC.h"
  36 #include "AliMathBase.h"
  37 #include "TTreeStream.h"
  38
  39 //date
  40 #include "event.h"
  41
  42 //header file
  43 #include "AliTPCCalibKr.h"
  44
  45 //----------------------------------------------------------------------------
  46 // The AliTPCCalibKr class description (TPC Kr calibration).
  47 //
  48 //
  49 // The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
  50 // its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).
  51 //
  52 // The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
  53 // These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration
  54 // parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
  55 // In addition the debugCalibKr.root file with debug information is created.
  56 //
  57
  58 /*
  59
  60 // Usage example:
  61 //
  62
  63 // 1. Analyse output histograms:
  64 TFile f("outHistFile.root");
  65 AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr");
  66 obj->SetRadius(0,0);
  67 obj->SetStep(1,1);
  68 obj->Analyse();
  69
  70 // 2. See calibration parameters e.g.:
  71 TFile f("calibKr.root");
  72 spectrMean->GetCalROC(70)->GetValue(40,40);
  73 fitMean->GetCalROC(70)->GetValue(40,40);
  74
  75 // 3. See debug information e.g.:
  76 TFile f("debugCalibKr.root");
  77 .ls;
  78
  79 // -- Print calibKr TTree content
  80 calibKr->Print();
  81
  82 // -- Draw calibKr TTree variables
  83 calibKr.Draw("fitMean");
  84
  85 */
  86
  87
  88 //
  89 // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
  90 //-----------------------------------------------------------------------------
  91
  92 ClassImp(AliTPCCalibKr)
  93
  94 AliTPCCalibKr::AliTPCCalibKr() :
  95   TObject(),
  96   fASide(kTRUE),
  97   fCSide(kTRUE),
  98   fHistoKrArray(72),
  99   fADCOverClustSizeMin(0.0),
 100   fADCOverClustSizeMax(1.0e9),
 101   fMaxADCOverClustADCMin(0.0),
 102   fMaxADCOverClustADCMax(1.0e9),
 103   fTimeMin(0.0),
 104   fTimeMax(1.0e9),
 105   fClustSizeMin(0.0),
 106   fClustSizeMax(1.0e9),
 107   fTimebinRmsIrocMin(0.0),
 108   fPadRmsIrocMin(0.0),
 109   fRowRmsIrocMin(0.0),
 110   fClusterPadSize1DIrocMax(200),
 111   fCurveCoefficientIroc(1.0e9),
 112   fTimebinRmsOrocMin(0.0),
 113   fPadRmsOrocMin(0.0),
 114   fRowRmsOrocMin(0.0),
 115   fClusterPadSize1DOrocMax(200),
 116   fCurveCoefficientOroc(1.0e9),
 117   fIrocHistogramMin(100.),
 118   fIrocHistogramMax(6000.),
 119   fIrocHistogramNbins(200),
 120   fOrocHistogramMin(100.),
 121   fOrocHistogramMax(5500.),
 122   fOrocHistogramNbins(200),
 123   fRowRadius(0),
 124   fPadRadius(0),
 125   fRowStep(1),
 126   fPadStep(1)
 127
 128 {
 129   //
 130   // default constructor
 131   //
 132
 133   // TObjArray with histograms
 134   fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
 135   fHistoKrArray.Clear();
 136
 137   // init cuts (by Stefan)
 138 //  SetADCOverClustSizeRange(7,200);
 139 //  SetMaxADCOverClustADCRange(0.01,0.4);
 140 //  SetTimeRange(200,600);
 141 //  SetClustSizeRange(6,200);
 142
 143   //init  cuts (by Adam)
 144   SetTimebinRmsMin(1.6,0.8);
 145   SetPadRmsMin(0.825,0.55);
 146   SetRowRmsMin(0.1,0.1);
 147   SetClusterPadSize1DMax(15,11);
 148   SetCurveCoefficient(1.,2.);
 149
 150   //set histograms settings
 151   SetIrocHistogram(200,100,6000);
 152   SetOrocHistogram(200,100,5500);
 153
 154   // init histograms
 155   //Init();
 156 }
 157
 158 //_____________________________________________________________________
 159 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
 160   TObject(pad),
 161
 162   fASide(pad.fASide),
 163   fCSide(pad.fCSide),
 164   fHistoKrArray(72),
 165   fADCOverClustSizeMin(pad.fADCOverClustSizeMin),
 166   fADCOverClustSizeMax(pad.fADCOverClustSizeMax),
 167   fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),
 168   fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),
 169   fTimeMin(pad.fTimeMin),
 170   fTimeMax(pad.fTimeMax),
 171   fClustSizeMin(pad.fClustSizeMin),
 172   fClustSizeMax(pad.fClustSizeMax),
 173   fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),
 174   fPadRmsIrocMin(pad.fPadRmsIrocMin),
 175   fRowRmsIrocMin(pad.fRowRmsIrocMin),
 176   fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),
 177   fCurveCoefficientIroc(pad.fCurveCoefficientIroc),
 178   fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),
 179   fPadRmsOrocMin(pad.fPadRmsOrocMin),
 180   fRowRmsOrocMin(pad.fRowRmsOrocMin),
 181   fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),
 182   fCurveCoefficientOroc(pad.fCurveCoefficientOroc),
 183   fIrocHistogramMin(pad.fIrocHistogramMin),
 184   fIrocHistogramMax(pad.fIrocHistogramMax),
 185   fIrocHistogramNbins(pad.fIrocHistogramNbins),
 186   fOrocHistogramMin(pad.fOrocHistogramMin),
 187   fOrocHistogramMax(pad.fOrocHistogramMax),
 188   fOrocHistogramNbins(pad.fOrocHistogramNbins),
 189   fRowRadius(pad.fRowRadius),
 190   fPadRadius(pad.fPadRadius),
 191   fRowStep(pad.fRowStep),
 192   fPadStep(pad.fPadStep)
 193
 194 {
 195   // copy constructor
 196
 197   // TObjArray with histograms
 198   fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
 199   fHistoKrArray.Clear();
 200
 201   for (Int_t iSec = 0; iSec < 72; ++iSec)
 202   {
 203     TH3F *hOld = pad.GetHistoKr(iSec);
 204         if(hOld) {
 205       TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
 206       fHistoKrArray.AddAt(hNew,iSec);
 207         }
 208   }
 209 }
 210
 211 //_____________________________________________________________________
 212 AliTPCCalibKr::~AliTPCCalibKr()
 213 {
 214   //
 215   // destructor
 216   //
 217
 218   // for (Int_t iSec = 0; iSec < 72; ++iSec) {
 219   //     if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
 220   // }
 221   fHistoKrArray.Delete();
 222 }
 223
 224 //_____________________________________________________________________
 225 AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)
 226 {
 227   // assignment operator
 228
 229   if (&source == this) return *this;
 230   new (this) AliTPCCalibKr(source);
 231
 232   return *this;
 233 }
 234
 235 //_____________________________________________________________________
 236 void AliTPCCalibKr::Init()
 237 {
 238   //
 239   // init output histograms
 240   //
 241
 242   // add histograms to the TObjArray
 243   for(Int_t i=0; i<72; ++i) {
 244
 245     // C - side
 246     if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
 247       TH3F *hist = CreateHisto(i);
 248       if(hist) fHistoKrArray.AddAt(hist,i);
 249     }
 250
 251     // A - side
 252     if(IsCSide(i) == kFALSE && fASide == kTRUE) {
 253       TH3F *hist = CreateHisto(i);
 254       if(hist) fHistoKrArray.AddAt(hist,i);
 255     }
 256   }
 257 }
 258
 259 //_____________________________________________________________________
 260 Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
 261 {
 262   //
 263   // process events
 264   // call event by event
 265   //
 266
 267   if(cluster) Update(cluster);
 268   else return kFALSE;
 269
 270  return kTRUE;
 271 }
 272
 273 //_____________________________________________________________________
 274 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
 275 {
 276     //
 277     // create new histogram
 278         //
 279     char name[256];
 280         TH3F *h;
 281
 282         snprintf(name,256,"ADCcluster_ch%d",chamber);
 283
 284     if( IsIROC(chamber) == kTRUE )
 285         {
 286           h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);
 287         } else {
 288           h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);
 289         }
 290     h->SetXTitle("padrow");
 291     h->SetYTitle("pad");
 292     h->SetZTitle("fADC");
 293
 294 return h;
 295 }
 296
 297 //_____________________________________________________________________
 298 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
 299 {
 300 // check if IROCs
 301 // returns kTRUE if IROCs and kFALSE if OROCs
 302
 303   if(chamber>=0 && chamber<36) return kTRUE;
 304
 305 return kFALSE;
 306 }
 307
 308 //_____________________________________________________________________
 309 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
 310 {
 311 // check if C side
 312 // returns kTRUE if C side and kFALSE if A side
 313
 314   if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
 315
 316 return kFALSE;
 317 }
 318
 319 //_____________________________________________________________________
 320 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)
 321 {
 322   //
 323   // fill existing histograms
 324   //
 325
 326   if (!Accept(cl)) return kFALSE;
 327   TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
 328   if(!h) return kFALSE;
 329
 330   h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
 331
 332   return kTRUE;
 333 }
 334
 335 //_____________________________________________________________________
 336 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){
 337   //
 338   // cuts
 339   //
 340   /*
 341     TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings -
 342     TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal
 343     TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal
 344     TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal
 345     TCut cutR4("cutR4","fMax.fTime>200");   // noise removal
 346     TCut cutR5("cutR5","fMax.fTime<600");   // noise removal
 347     TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks
 348     TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;
 349   */
 350   /*
 351   //R0
 352   if ((float)cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;
 353   // R1
 354   if ((float)cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;
 355   //R2
 356   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;
 357   //R3
 358   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
 359   //R4
 360   if (cl->GetMax().GetTime() < 200) return kFALSE;
 361   //R5
 362   if (cl->GetMax().GetTime() > 600) return kFALSE;
 363   //S1
 364   if (cl->GetSize()>200) return kFALSE;
 365   if (cl->GetSize()<6)  return kFALSE;
 366
 367   SetADCOverClustSizeRange(7,200);
 368   SetMaxADCOverClustADCRange(0.01,0.4);
 369   SetTimeRange(200,600);
 370   SetClustSizeRange(6,200);
 371   */
 372
 373   //R0
 374   if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax)        return kFALSE;
 375   // R1
 376   if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin)          return kFALSE;
 377   //R2
 378   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax)  return kFALSE;
 379   //R3
 380   if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
 381   //R4
 382   if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
 383   //R5
 384   if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
 385   //S1
 386   if (cl->GetSize()>fClustSizeMax) return kFALSE;
 387   if (cl->GetSize()<fClustSizeMin)  return kFALSE;
 388
 389   //
 390   // cuts by Adam
 391   //
 392   /*
 393     TCut cutAI0("cutAI0","fTimebinRMS>1.6");
 394     TCut cutAI1("cutAI1","fPadRMS>0.825");
 395     TCut cutAI2("cutAI2","fRowRMS>0.1");
 396     TCut cutAI3("cutAI3","fPads1D<15");
 397     TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0");
 398
 399     TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;
 400
 401     TCut cutAO0("cutAO0","fTimebinRMS>0.8");
 402     TCut cutAO1("cutAO1","fPadRMS>0.55");
 403     TCut cutAO2("cutAO2","fRowRMS>0.1");
 404     TCut cutAO3("cutAO3","fPads1D<11");
 405     TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0");
 406
 407     TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;
 408     TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");
 409
 410   */
 411   if(cl->GetSec()<36){  //IROCs
 412     //AI0
 413     if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;
 414     //AI1
 415     if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;
 416     //AI2
 417     if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;
 418     //AI3
 419     if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;
 420     //AI4
 421     if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
 422   }else{  //OROCs
 423     //AO0
 424     if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;
 425     //AO1
 426     if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;
 427     //AO2
 428     if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;
 429     //AO3
 430     if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;
 431     //AO4
 432     if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
 433   }
 434
 435   return kTRUE;
 436
 437 }
 438
 439 //_____________________________________________________________________
 440 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
 441 {
 442   // get histograms from fHistoKrArray
 443   return (TH3F*) fHistoKrArray.At(chamber);
 444 }
 445
 446 //_____________________________________________________________________
 447 void AliTPCCalibKr::Analyse()
 448 {
 449   //
 450   // analyse the histograms and extract krypton calibration parameters
 451   //
 452
 453   // AliTPCCalPads that will contain the calibration parameters
 454   AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
 455   AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
 456   AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
 457   AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
 458   AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
 459   AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
 460
 461   // file stream for debugging purposes
 462   TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
 463
 464   // if entries in spectrum less than minEntries, then the fit won't be performed
 465   Int_t minEntries = 1; //300;
 466
 467   Double_t windowFrac = 0.12;
 468   // the 3d histogram will be projected on the pads given by the following window size
 469   // set the numbers to 0 if you want to do a pad-by-pad calibration
 470   UInt_t rowRadius = fRowRadius;//4;
 471   UInt_t padRadius = fPadRadius;//4;
 472
 473   // the step size by which pad and row are incremented is given by the following numbers
 474   // set them to 1 if you want the finest granularity
 475   UInt_t rowStep = fRowStep;//1;//2    // formerly: 2*rowRadius
 476   UInt_t padStep = fPadStep;//1;//4    // formerly: 2*padRadius
 477
 478   for (Int_t chamber = 0; chamber < 72; chamber++) {
 479     //if (chamber != 71) continue;
 480     AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()
 481
 482     // Usually I would traverse each pad, take the spectrum for its neighbourhood and
 483     // obtain the calibration parameters. This takes very long, so instead I assign the same
 484     // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
 485     UInt_t nRows = roc.GetNrows();
 486     for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
 487       UInt_t nPads = roc.GetNPads(iRow);
 488       //if (iRow >= 10) break;
 489       for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
 490         //if (iPad >= 20) break;
 491         TH3F* h = GetHistoKr(chamber);
 492         if (!h) continue;
 493
 494         // the 3d histogram will be projected on the pads given by the following bounds
 495         // for rows and pads
 496         Int_t rowLow = iRow - rowRadius;
 497         UInt_t rowUp = iRow + rowRadius + rowStep-1;
 498         Int_t padLow = iPad - padRadius;
 499         UInt_t padUp = iPad + padRadius + padStep-1;
 500         // if window goes out of chamber
 501         if (rowLow < 0) rowLow = 0;
 502         if (rowUp >= nRows) rowUp = nRows - 1;
 503         if (padLow < 0) padLow = 0;
 504         if (padUp >= nPads) padUp = nPads - 1;
 505
 506         // project the histogram
 507         //TH1D* projH = h->ProjectionZ("projH", rowLow+1, rowUp+1, padLow+1, padUp+1); // SLOW
 508         TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);
 509
 510         // get the number of entries in the spectrum
 511         Double_t entries = projH->GetEntries();
 512         if (entries < minEntries) { delete projH; continue; }
 513
 514         // get the two calibration parameters mean of spectrum and RMS of spectrum
 515         Double_t histMean = projH->GetMean();
 516         Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
 517
 518         // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
 519         Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
 520                 Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);
 521                 Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
 522         Double_t integCharge =  projH->Integral(minBin,maxBin);
 523
 524         Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
 525
 526         if (fitResult != 0) {
 527           Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
 528           //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f,  %f\n", chamber, iRow, iPad, entries, maxEntries);
 529
 530           delete projH;
 531           continue;
 532         }
 533
 534         // get the two calibration parameters mean of gauss fit and sigma of gauss fit
 535         TF1* gausFit = projH->GetFunction("gaus");
 536         Double_t fitMean = gausFit->GetParameter(1);
 537         Double_t fitRMS = gausFit->GetParameter(2);
 538         Int_t numberFitPoints = gausFit->GetNumberFitPoints();
 539         if (numberFitPoints == 0) continue;
 540         Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
 541         delete gausFit;
 542         delete projH;
 543         if (fitMean <= 0) continue;
 544         //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
 545
 546         // write the calibration parameters for each pad that the 3d histogram was projected onto
 547         // (with considering the step size) to the CalPads
 548         // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
 549         // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
 550         for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
 551           if (r < 0 || r >= (Int_t)nRows) continue;
 552           UInt_t nPadsR = roc.GetNPads(r);
 553           for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
 554             if (p < 0 || p >= (Int_t)nPadsR) continue;
 555             spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
 556             spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
 557             fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
 558             fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
 559             fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
 560             entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
 561
 562             (*debugStream) << "calibKr" <<
 563               "sector=" << chamber <<          // chamber number
 564               "row=" << r <<                   // row number
 565               "pad=" << p <<                   // pad number
 566               "histMean=" << histMean <<       // mean of the spectrum
 567               "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean
 568               "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak
 569               "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak
 570               "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit
 571               "entries=" << entries <<         // number of entries for the spectrum
 572               "\n";
 573           }
 574         }
 575       }
 576     }
 577   }
 578
 579   TFile f("calibKr.root", "recreate");
 580   spectrMeanCalPad->Write();
 581   spectrRMSCalPad->Write();
 582   fitMeanCalPad->Write();
 583   fitRMSCalPad->Write();
 584   fitNormChi2CalPad->Write();
 585   entriesCalPad->Write();
 586   f.Close();
 587   delete spectrMeanCalPad;
 588   delete spectrRMSCalPad;
 589   delete fitMeanCalPad;
 590   delete fitRMSCalPad;
 591   delete fitNormChi2CalPad;
 592   delete entriesCalPad;
 593   delete debugStream;
 594 }
 595
 596 //_____________________________________________________________________
 597 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
 598 {
 599   // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
 600   // replaces TH3F::ProjectZ() to gain more speed
 601
 602   TAxis* xAxis = histo3D->GetXaxis();
 603   TAxis* yAxis = histo3D->GetYaxis();
 604   TAxis* zAxis = histo3D->GetZaxis();
 605   Double_t zMinVal = zAxis->GetXmin();
 606   Double_t zMaxVal = zAxis->GetXmax();
 607
 608   Int_t nBinsZ = zAxis->GetNbins();
 609   TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
 610
 611   Int_t nx = xAxis->GetNbins()+2;
 612   Int_t ny = yAxis->GetNbins()+2;
 613   Int_t bin = 0;
 614   Double_t entries = 0.;
 615   for (Int_t x = xMin; x <= xMax; x++) {
 616     for (Int_t y = yMin; y <= yMax; y++) {
 617       for (Int_t z = 0; z <= nBinsZ+1; z++) {
 618         bin = x + nx * (y + ny * z);
 619         Double_t val = histo3D->GetBinContent(bin);
 620         projH->Fill(zAxis->GetBinCenter(z), val);
 621         entries += val;
 622       }
 623     }
 624   }
 625   projH->SetEntries((Long64_t)entries);
 626   return projH;
 627 }
 628
 629 //_____________________________________________________________________
 630 Long64_t AliTPCCalibKr::Merge(TCollection* list) {
 631 // merge function
 632 //
 633
 634 if (!list)
 635 return 0;
 636
 637 if (list->IsEmpty())
 638 return 1;
 639
 640 TIterator* iter = list->MakeIterator();
 641 TObject* obj = 0;
 642
 643     iter->Reset();
 644     Int_t count=0;
 645         while((obj = iter->Next()) != 0)
 646         {
 647           AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
 648           if (entry == 0) continue;
 649
 650                 for(int i=0; i<72; ++i) {
 651                   if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
 652                   ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
 653                   }
 654
 655                   if(IsCSide(i) == kFALSE && fASide == kTRUE) {
 656                     ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
 657                   }
 658                 }
 659
 660           count++;
 661         }
 662
 663 return count;
 664 }