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