TPC/AliTPCCalibKr.cxx

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