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