TPC/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 "AliTPCRawStream.h"
  34 #include "AliTPCCalROC.h"
  35 #include "AliTPCCalPad.h"
  36 #include "AliTPCROC.h"
  37 #include "AliMathBase.h"
  38 #include "TTreeStream.h"
  39 #include "AliTPCRawStreamFast.h"
  40
  41 //date
  42 #include "event.h"
  43
  44 //header file
  45 #include "AliTPCCalibKr.h"
  46
  47 //----------------------------------------------------------------------------
  48 // The AliTPCCalibKr class description (TPC Kr calibration).
  49 //
  50 //
  51 // The AliTPCCalibKr fills the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
  52 // its data memebers.
  53 //
  54 // As the input it requires the tree with reconstructed Kr clusters (AliTPCclusterKr objects).
  55 // The AliTPCCalibKr objects containing an array of TH3F histograms are stored (by default) in the
  56 // ouptut (outHistFile.root) file.
  57 //
  58 // The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
  59 // These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration
  60 // parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
  61 // In addition the debugCalibKr.root file with debug information is created.
  62 //
  63 // Usage example:
  64 //
  65 // 1. Create outHistFile.root histogram file:
  66 //
  67 // -- Load libXrdClient.so if data on Xrd cluster e.g. (GSI)
  68 // gSystem->Load("/usr/local/grid/XRootd/GSI/lib64/libXrdClient.so");
  69 //
  70 // -- Load toolkit
  71 // gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  72 // gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
  73 // AliXRDPROOFtoolkit tool;
  74 //
  75 // -- Make chain of files
  76 // TChain * chain = tool.MakeChain("KrClusters.txt","Kr","",1000,0);
  77 //
  78 // -- Run AliTPCCalibKr task (Only TPC C side)
  79 // AliTPCCalibKr *task = new AliTPCCalibKr;
  80 // task->SetInputChain(chain);
  81 // task->SetASide(kFALSE);
  82 //
  83 // task->Process();
  84 //
  85 // 2. Analyse output histograms:
  86 //
  87 // TFile f("outHistFile.root");
  88 // AliTPCCalibKr.Analyse();
  89 //
  90 // 3. See calibration parameters e.g.:
  91 //
  92 // TFile f("calibKr.root");
  93 // spectrMean->GetCalROC(70)->GetValue(40,40);
  94 // fitMean->GetCalROC(70)->GetValue(40,40);
  95 //
  96 // 4. See debug information e.g.:
  97 //
  98 // TFile f("debugCalibKr.root");
  99 // .ls;
 100 //
 101 // -- Print calibKr TTree content
 102 // calibKr->Print();
 103 //
 104 // -- Draw calibKr TTree variables
 105 // calibKr.Draw("fitMean");
 106 //
 107 //
 108 // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
 109 //-----------------------------------------------------------------------------
 110
 111 ClassImp(AliTPCCalibKr)
 112
 113 AliTPCCalibKr::AliTPCCalibKr() :
 114   TObject(),
 115
 116   bOutputHisto(kTRUE),
 117   bASide(kTRUE),
 118   bCSide(kTRUE),
 119   fClusters(0),
 120   fClustKr(0),
 121   fTree(0),
 122   fHistoKrArray(72)
 123 {
 124   //
 125   // default constructor
 126   //
 127 }
 128
 129 //_____________________________________________________________________
 130 AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
 131   TObject(pad),
 132
 133   bOutputHisto(pad.bOutputHisto),
 134   bASide(pad.bASide),
 135   bCSide(pad.bCSide),
 136   fClusters(pad.fClusters),
 137   fClustKr(pad.fClustKr),
 138   fTree(pad.fTree),
 139   fHistoKrArray(72)
 140 {
 141   // copy constructor
 142
 143   for (Int_t iSec = 0; iSec < 72; ++iSec)
 144   {
 145     TH3F *hOld = pad.GetHistoKr(iSec);
 146         if(hOld) {
 147       TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
 148       fHistoKrArray.AddAt(hNew,iSec);
 149         }
 150   }
 151 }
 152
 153 //_____________________________________________________________________
 154 AliTPCCalibKr::~AliTPCCalibKr()
 155 {
 156   //
 157   // destructor
 158   //
 159   if(fClustKr)  delete fClustKr; fClustKr = 0;
 160   if(fClusters) delete fClusters; fClusters = 0;
 161   if(fTree)     delete fTree; fTree = 0;
 162   fHistoKrArray.Delete();
 163 }
 164
 165 //_____________________________________________________________________
 166 AliTPCCalibKr& AliTPCCalibKr::operator = (const  AliTPCCalibKr &source)
 167 {
 168   // assignment operator
 169
 170   if (&source == this) return *this;
 171   new (this) AliTPCCalibKr(source);
 172
 173   return *this;
 174 }
 175
 176 //_____________________________________________________________________
 177 void AliTPCCalibKr::Init()
 178 {
 179   //
 180   // init input tree and output histograms
 181   //
 182
 183   // set input tree
 184   if(!fTree) {
 185    Printf("ERROR: Could not read chain from input");
 186   }
 187   else {
 188    fTree->SetBranchStatus("*",1);
 189   }
 190
 191   // set branch address
 192   fClusters = new TClonesArray("AliTPCclusterKr");
 193
 194   if(!fTree->GetBranch("fClusters")) {
 195     Printf("ERROR: Could not get fClusters branch from input");
 196   } else {
 197    fTree->GetBranch("fClusters")->SetAddress(&fClusters);
 198   }
 199
 200   // create output TObjArray
 201   fHistoKrArray.Clear();
 202
 203   // add histograms to the TObjArray
 204   for(Int_t i=0; i<72; ++i) {
 205
 206         // C - side
 207         if( IsCSide(i) == kTRUE && bCSide == kTRUE) {
 208       TH3F *hist = CreateHisto(i);
 209       if(hist) fHistoKrArray.AddAt(hist,i);
 210         }
 211
 212         // A - side
 213         if(IsCSide(i) == kFALSE && bASide == kTRUE) {
 214       TH3F *hist = CreateHisto(i);
 215       if(hist) fHistoKrArray.AddAt(hist,i);
 216         }
 217
 218   }
 219 }
 220
 221 //_____________________________________________________________________
 222 Bool_t AliTPCCalibKr::ReadEntry(Int_t evt)
 223 {
 224   //
 225   // read entry from the tree
 226   //
 227   Long64_t centry = fTree->LoadTree(evt);
 228   if(centry < 0) return kFALSE;
 229
 230   if(!fTree->GetBranch("fClusters"))
 231   {
 232     Printf("ERROR: Could not get fClusters branch from input");
 233         return kFALSE;
 234   } else {
 235    fTree->GetBranch("fClusters")->SetAddress(&fClusters);
 236   }
 237
 238   fTree->GetEntry(evt);
 239
 240 return kTRUE;
 241 }
 242
 243 //_____________________________________________________________________
 244 Bool_t AliTPCCalibKr::Process()
 245 {
 246   //
 247   // process events
 248   // call event by event
 249   //
 250
 251   // init tree
 252   Init();
 253
 254   // get events
 255   if(!fTree) return kFALSE;
 256   Int_t nEvents = fTree->GetEntries();
 257
 258   // fill histograms
 259   for(Int_t i=0; i<nEvents; ++i)
 260   {
 261     if(ReadEntry(i) == kFALSE) return kFALSE;
 262
 263     if(!(i%10000)) cout << "evt: " << i << endl;
 264
 265     // get TClonesArray entries
 266     fClustKr = 0;
 267     Int_t entries = fClusters->GetEntries();
 268     for(Int_t j=0; j < entries; ++j)
 269         {
 270           fClustKr = (AliTPCclusterKr*)fClusters->At(j);
 271
 272       if(fClustKr) Update(fClustKr);
 273           else return kFALSE;
 274         }
 275   }
 276
 277   // write output
 278   return Terminate();
 279 }
 280
 281 //_____________________________________________________________________
 282 TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
 283 {
 284     //
 285     // create new histogram
 286         //
 287     char name[256];
 288         TH3F *h;
 289
 290     sprintf(name,"ADCcluster_ch%d",chamber);
 291
 292     if( IsIROC(chamber) == kTRUE )
 293         {
 294            h = new TH3F(name,name,63,0,63,100,0,100,150,100,3000);
 295         } else {
 296            h = new TH3F(name,name,96,0,96,100,0,100,150,100,3000);
 297         }
 298     h->SetXTitle("padrow");
 299     h->SetYTitle("pad");
 300     h->SetZTitle("fADC");
 301
 302 return h;
 303 }
 304
 305 //_____________________________________________________________________
 306 Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
 307 {
 308 // check if IROCs
 309 // returns kTRUE if IROCs and kFALSE if OROCs
 310
 311   if(chamber>=0 && chamber<36) return kTRUE;
 312
 313 return kFALSE;
 314 }
 315
 316 //_____________________________________________________________________
 317 Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
 318 {
 319 // check if C side
 320 // returns kTRUE if C side and kFALSE if A side
 321
 322   if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
 323
 324 return kFALSE;
 325 }
 326
 327 //_____________________________________________________________________
 328 Bool_t AliTPCCalibKr::Update(AliTPCclusterKr  *cl)
 329 {
 330   //
 331   // fill existing histograms
 332   //
 333
 334   if (!Accept(cl)) return kFALSE;
 335   TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
 336   if(!h) return kFALSE;
 337
 338   h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
 339
 340   return kTRUE;
 341 }
 342
 343 Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr  *cl){
 344   //
 345   // cuts
 346   //
 347   /*
 348     TCut cutR0("cutR0","fADCcluster/fSize<200");        // adjust it according v seetings -
 349     TCut cutR1("cutR1","fADCcluster/fSize>7");          // cosmic tracks and noise removal
 350     TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4");    // digital noise removal
 351     TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01");   // noise removal
 352     TCut cutS1("cutS1","fSize<200");    // adjust it according v seetings - cosmic tracks
 353     TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
 354   */
 355   //R0
 356   if (cl->GetADCcluster()/ cl->GetSize() >200)        return kFALSE;
 357   // R1
 358   if (cl->GetADCcluster()/ cl->GetSize() <7)          return kFALSE;
 359   //R2
 360   if (cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4)  return kFALSE;
 361   //R3
 362   if (cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
 363   //S1
 364   if (cl->GetSize()>200) return kFALSE;
 365   if (cl->GetSize()<6)  return kFALSE;
 366   return kTRUE;
 367
 368 }
 369
 370
 371
 372 //_____________________________________________________________________
 373 TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
 374 {
 375   // get histograms from fHistoKrArray
 376   return (TH3F*) fHistoKrArray.At(chamber);
 377 }
 378
 379 //_____________________________________________________________________
 380 Bool_t AliTPCCalibKr::Terminate()
 381 {
 382   //
 383   // store AliTPCCalibKr in the output file
 384   //
 385   if(bOutputHisto) {
 386     TFile *outFile = new TFile("outHistFile.root","RECREATE");
 387
 388     if(outFile)
 389         {
 390           outFile->cd();
 391
 392           for(int i=0; i<72; ++i) {
 393              if( IsCSide(i) == kTRUE && bCSide == kTRUE)
 394                printf("C side chamber: %d, 3D histo entries: %10.f \n",i,((TH3F*)fHistoKrArray.At(i))->GetEntries());
 395
 396              if( IsCSide(i) == kFALSE && bASide == kTRUE)
 397                printf("A side chamber: %d, 3D histo entries: %10.f \n",i,((TH3F*)fHistoKrArray.At(i))->GetEntries());
 398           }
 399           this->Write();
 400           outFile->Close();
 401
 402           return kTRUE;
 403         }
 404         else
 405           return kFALSE;
 406   }
 407
 408 return kFALSE;
 409 }
 410
 411 //_____________________________________________________________________
 412 void AliTPCCalibKr::Analyse()
 413 {
 414   //
 415   // analyse the histograms and extract krypton calibration parameters
 416   //
 417
 418   // AliTPCCalPads that will contain the calibration parameters
 419   AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
 420   AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
 421   AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
 422   AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
 423   AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
 424   AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
 425
 426   // file stream for debugging purposes
 427   TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
 428
 429   // if entries in spectrum less than minEntries, then the fit won't be performed
 430   Int_t minEntries = 1; //300;
 431
 432   Double_t windowFrac = 0.12;
 433   // the 3d histogram will be projected on the pads given by the following window size
 434   // set the numbers to 0 if you want to do a pad-by-pad calibration
 435   UInt_t rowRadius = 5;
 436   UInt_t padRadius = 10;
 437   // the step size by which pad and row are incremented is given by the following numbers
 438   // set them to 1 if you want the finest granularity
 439   UInt_t rowStep = 1;    // formerly: 2*rowRadius
 440   UInt_t padStep = 1;    // formerly: 2*padRadius
 441
 442   for (Int_t chamber = 0; chamber < 72; chamber++) {
 443     //if (chamber != 71) continue;
 444     AliTPCCalROC roc(chamber);  // I need this only for GetNrows() and GetNPads()
 445
 446     // Usually I would traverse each pad, take the spectrum for its neighbourhood and
 447     // obtain the calibration parameters. This takes very long, so instead I assign the same
 448     // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
 449     UInt_t nRows = roc.GetNrows();
 450     for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
 451       UInt_t nPads = roc.GetNPads(iRow);
 452       //if (iRow >= 10) break;
 453       for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
 454         //if (iPad >= 20) break;
 455         TH3F* h = GetHistoKr(chamber);
 456         if (!h) continue;
 457
 458         // the 3d histogram will be projected on the pads given by the following bounds
 459         // for rows and pads
 460         Int_t rowLow = iRow - rowRadius;
 461         UInt_t rowUp = iRow + rowRadius;
 462         Int_t padLow = iPad - padRadius;
 463         UInt_t padUp = iPad + padRadius;
 464         // if window goes out of chamber
 465         if (rowLow < 0) rowLow = 0;
 466         if (rowUp >= nRows) rowUp = nRows - 1;
 467         if (padLow < 0) padLow = 0;
 468         if (padUp >= nPads) padUp = nPads - 1;
 469
 470         // project the histogram
 471         //TH1D* projH = h->ProjectionZ("projH", rowLow, rowUp, padLow, padUp); // SLOW
 472         TH1D* projH = ProjectHisto(h, "projH", rowLow, rowUp, padLow, padUp);
 473
 474         // get the number of entries in the spectrum
 475         Double_t entries = projH->GetEntries();
 476         if (entries < minEntries) { delete projH; continue; }
 477
 478         // get the two calibration parameters mean of spectrum and RMS of spectrum
 479         Double_t histMean = projH->GetMean();
 480         Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
 481
 482         // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
 483         Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
 484         Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
 485         if (fitResult != 0) {
 486           Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i.", chamber, rowLow, rowUp, padLow, padUp);
 487           delete projH;
 488           continue;
 489         }
 490
 491         // get the two calibration parameters mean of gauss fit and sigma of gauss fit
 492         TF1* gausFit = projH->GetFunction("gaus");
 493         Double_t fitMean = gausFit->GetParameter(1);
 494         Double_t fitRMS = gausFit->GetParameter(2);
 495         Int_t numberFitPoints = gausFit->GetNumberFitPoints();
 496         if (numberFitPoints == 0) continue;
 497         Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
 498         delete gausFit;
 499         delete projH;
 500         if (fitMean <= 0) continue;
 501         printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
 502
 503         // write the calibration parameters for each pad that the 3d histogram was projected onto
 504         // (with considering the step size) to the CalPads
 505         // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
 506         // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
 507         for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
 508           if (r < 0 || r >= (Int_t)nRows) continue;
 509           UInt_t nPads = roc.GetNPads(r);
 510           for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
 511             if (p < 0 || p >= (Int_t)nPads) continue;
 512             spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
 513             spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
 514             fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
 515             fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
 516             fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
 517             entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
 518
 519             (*debugStream) << "calibKr" <<
 520               "sector=" << chamber <<          // chamber number
 521               "row=" << r <<                   // row number
 522               "pad=" << p <<                   // pad number
 523               "histMean=" << histMean <<       // mean of the spectrum
 524               "histRMS=" << histRMS <<         // RMS of the spectrum divided by the mean
 525               "fitMean=" << fitMean <<         // Gauss fitted mean of the 41.6 keV Kr peak
 526               "fitRMS=" << fitRMS <<           // Gauss fitted sigma of the 41.6 keV Kr peak
 527               "fitNormChi2" << fitNormChi2 <<  // normalized chi square of the Gauss fit
 528               "entries=" << entries <<         // number of entries for the spectrum
 529               "\n";
 530           }
 531         }
 532       }
 533     }
 534   }
 535
 536   TFile f("calibKr.root", "recreate");
 537   spectrMeanCalPad->Write();
 538   spectrRMSCalPad->Write();
 539   fitMeanCalPad->Write();
 540   fitRMSCalPad->Write();
 541   fitNormChi2CalPad->Write();
 542   entriesCalPad->Write();
 543   f.Close();
 544   delete spectrMeanCalPad;
 545   delete spectrRMSCalPad;
 546   delete fitMeanCalPad;
 547   delete fitRMSCalPad;
 548   delete fitNormChi2CalPad;
 549   delete entriesCalPad;
 550   delete debugStream;
 551 }
 552
 553 //_____________________________________________________________________
 554 TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
 555 {
 556   // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
 557   // replaces TH3F::ProjectZ() to gain more speed
 558
 559   TAxis* xAxis = histo3D->GetXaxis();
 560   TAxis* yAxis = histo3D->GetYaxis();
 561   TAxis* zAxis = histo3D->GetZaxis();
 562   Double_t zMinVal = zAxis->GetXmin();
 563   Double_t zMaxVal = zAxis->GetXmax();
 564
 565   Int_t nBinsZ = zAxis->GetNbins();
 566   TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
 567
 568   Int_t nx = xAxis->GetNbins()+2;
 569   Int_t ny = yAxis->GetNbins()+2;
 570   Int_t bin = 0;
 571   Double_t entries = 0.;
 572   for (Int_t x = xMin; x <= xMax; x++) {
 573     for (Int_t y = yMin; y <= yMax; y++) {
 574       for (Int_t z = 0; z <= nBinsZ+1; z++) {
 575         bin = x + nx * (y + ny * z);
 576         Double_t val = histo3D->GetBinContent(bin);
 577         projH->Fill(zAxis->GetBinCenter(z), val);
 578         entries += val;
 579       }
 580     }
 581   }
 582   projH->SetEntries((Long64_t)entries);
 583   return projH;
 584 }