TPC/AliTPCcalibBase.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 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  Base class for the calibration components using
  20 //  as input TPCseeds and ESDs
  21 //  Event loop outside of the component
  22 //
  23 //
  24 // Base functionality to be implemeneted by component
  25 /*
  26    //In some cases only one of this function to be implemented
  27    virtual void     Process(AliESDEvent *event)
  28    virtual void     Process(AliTPCseed *track)
  29    //
  30    virtual Long64_t Merge(TCollection *li);
  31    virtual void     Analyze()
  32    void             Terminate();
  33 */
  34 // Functionality provided by base class for Algorith debuging:
  35 //  TTreeSRedirector * cstream =  GetDebugStreamer() - get debug streamer which can be use for numerical debugging
  36 //
  37
  38
  39
  40 //  marian.ivanov@cern.ch
  41 //
  42 #include "AliTPCcalibBase.h"
  43 #include "TSystem.h"
  44 #include "TFile.h"
  45 #include "TTreeStream.h"
  46 #include "TTimeStamp.h"
  47 #include "TGraph.h"
  48 #include "TGraphErrors.h"
  49 #include "TF1.h"
  50 #include "TH1.h"
  51 #include "THnSparse.h"
  52 #include "TH1D.h"
  53 #include "TH2D.h"
  54 #include "TAxis.h"
  55 #include "AliRecoParam.h"
  56
  57
  58 #include "AliLog.h"
  59 #include "AliESDEvent.h"
  60
  61
  62 ClassImp(AliTPCcalibBase)
  63
  64 AliTPCcalibBase::AliTPCcalibBase():
  65     TNamed(),
  66     fDebugStreamer(0),
  67     fStreamLevel(0),
  68     fRun(0),                  //!  current Run number
  69     fEvent(0),                //!  current Event number
  70     fTime(0),                 //!  current Time
  71     fTrigger(0),              //! current trigger type
  72     fMagF(0),                 //! current magnetic field
  73     fTriggerMaskReject(-1),   //trigger mask - reject trigger
  74     fTriggerMaskAccept(-1),   //trigger mask - accept trigger
  75     fHasLaser(kFALSE),                    //flag the laser is overlayed with given event
  76     fRejectLaser(kTRUE),                 //flag- reject laser
  77     fTriggerClass(),
  78     fCurrentEvent(0),           //! current event
  79     fCurrentTrack(0),           //! current esd track
  80     fCurrentFriendTrack(0),           //! current esd track
  81     fCurrentSeed(0),            //! current seed
  82     fDebugLevel(0)
  83 {
  84   //
  85   // Constructor
  86   //
  87 }
  88
  89 AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
  90   TNamed(name,title),
  91   fDebugStreamer(0),
  92   fStreamLevel(0),
  93   fRun(0),                  //!  current Run number
  94   fEvent(0),                //!  current Event number
  95   fTime(0),                 //!  current Time
  96   fTrigger(0),              //! current trigger type
  97   fMagF(0),                 //! current magnetic field
  98   fTriggerMaskReject(-1),   //trigger mask - reject trigger
  99   fTriggerMaskAccept(-1),   //trigger mask - accept trigger
 100   fHasLaser(kFALSE),                    //flag the laser is overlayed with given event
 101   fRejectLaser(kTRUE),                 //flag- reject laser
 102   fTriggerClass(),
 103   fCurrentEvent(0),           //! current event
 104   fCurrentTrack(0),           //! current esd track
 105   fCurrentFriendTrack(0),           //! current esd track
 106   fCurrentSeed(0),            //! current seed
 107   fDebugLevel(0)
 108 {
 109   //
 110   // Constructor
 111   //
 112 }
 113
 114 AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
 115   TNamed(calib),
 116   fDebugStreamer(0),
 117   fStreamLevel(calib.fStreamLevel),
 118   fRun(0),                  //!  current Run number
 119   fEvent(0),                //!  current Event number
 120   fTime(0),                 //!  current Time
 121   fTrigger(0),              //! current trigger type
 122   fMagF(0),                 //! current magnetic field
 123   fTriggerMaskReject(calib.fTriggerMaskReject),   //trigger mask - reject trigger
 124   fTriggerMaskAccept(calib.fTriggerMaskAccept),   //trigger mask - accept trigger
 125   fHasLaser(calib.fHasLaser),                    //flag the laser is overlayed with given event
 126   fRejectLaser(calib.fRejectLaser),                 //flag- reject laser
 127   fTriggerClass(calib.fTriggerClass),
 128   fCurrentEvent(0),           //! current event
 129   fCurrentTrack(0),           //! current esd track
 130   fCurrentFriendTrack(0),           //! current esd track
 131   fCurrentSeed(0),            //! current seed
 132   fDebugLevel(calib.fDebugLevel)
 133 {
 134   //
 135   // copy constructor
 136   //
 137 }
 138
 139 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
 140   //
 141   // operator=
 142   //
 143   if (this == &calib) return (*this);
 144
 145   ((TNamed *)this)->operator=(calib);
 146   fDebugStreamer=0;
 147   fStreamLevel=calib.fStreamLevel;
 148   fDebugLevel=calib.fDebugLevel;
 149   return *this;
 150 }
 151
 152
 153 AliTPCcalibBase::~AliTPCcalibBase() {
 154   //
 155   // destructor
 156   //
 157   if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
 158   if (fDebugStreamer) delete fDebugStreamer;
 159   fDebugStreamer=0;
 160 }
 161
 162 void  AliTPCcalibBase::Terminate(){
 163   //
 164   //
 165   //
 166   if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
 167   if (fDebugStreamer) delete fDebugStreamer;
 168   fDebugStreamer = 0;
 169   return;
 170 }
 171
 172 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
 173   //
 174   // Get Debug streamer
 175   // In case debug streamer not yet initialized and StreamLevel>0 create new one
 176   //
 177   if (fStreamLevel==0) return 0;
 178   if (fDebugStreamer) return fDebugStreamer;
 179   TString dsName;
 180   dsName=GetName();
 181   dsName+="Debug.root";
 182   dsName.ReplaceAll(" ","");
 183   fDebugStreamer = new TTreeSRedirector(dsName.Data());
 184   return fDebugStreamer;
 185 }
 186
 187
 188 void    AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
 189   //
 190   //
 191   //
 192   fRun     = event->GetRunNumber();
 193   fEvent   = event->GetEventNumberInFile();
 194   fTime    = event->GetTimeStamp();
 195   fTrigger = event->GetTriggerMask();
 196   fMagF    = event->GetMagneticField();
 197   fTriggerClass = event->GetFiredTriggerClasses().Data();
 198   fHasLaser = HasLaser(event);
 199
 200 }
 201
 202
 203 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
 204   //
 205   //
 206   //
 207   // Use event specie
 208   Bool_t isLaser=kFALSE;
 209   UInt_t specie = event->GetEventSpecie();  // select only cosmic events
 210   if (specie==AliRecoParam::kCalib) {
 211     isLaser = kTRUE;
 212   }
 213   return isLaser;
 214 }
 215
 216
 217
 218 Bool_t AliTPCcalibBase::AcceptTrigger(){
 219   //
 220   // Apply trigger mask - Don't do calibration for non proper triggers
 221   //
 222   if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
 223   if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
 224   if (fHasLaser && fRejectLaser) return kFALSE;
 225   return kTRUE;
 226 }
 227
 228
 229 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
 230   //
 231   // store  - copy debug output to the destination position
 232   // currently ONLY for local copy
 233   if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
 234   if (fStreamLevel==0) return;
 235   TString dsName;
 236   dsName=GetName();
 237   dsName+="Debug.root";
 238   dsName.ReplaceAll(" ","");
 239   TString dsName2=path;
 240   gSystem->MakeDirectory(dsName2.Data());
 241   dsName2+=gSystem->HostName();
 242   gSystem->MakeDirectory(dsName2.Data());
 243   dsName2+="/";
 244   TTimeStamp s;
 245   dsName2+=Int_t(s.GetNanoSec());
 246   dsName2+="/";
 247   gSystem->MakeDirectory(dsName2.Data());
 248   dsName2+=dsName;
 249   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
 250   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
 251   TFile::Cp(dsName.Data(),dsName2.Data());
 252 }
 253
 254
 255
 256 TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp, Bool_t useMedian, TTreeSRedirector *cstream, Int_t ival){
 257   //
 258   // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
 259   //
 260   TF1 funcGaus("funcGaus","gaus");
 261   TH2D * hist = h->Projection(axisDim1, axisDim2);
 262   Double_t *xvec = new Double_t[hist->GetNbinsX()];
 263   Double_t *yvec = new Double_t[hist->GetNbinsX()];
 264   Double_t *xerr = new Double_t[hist->GetNbinsX()];
 265   Double_t *yerr = new Double_t[hist->GetNbinsX()];
 266   Int_t counter  = 0;
 267   TH1D * projectionHist =0;
 268   //
 269
 270   for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
 271     Int_t nsum=0;
 272     Int_t imin   =  i;
 273     Int_t imax   =  i;
 274
 275     for (Int_t idelta=0; idelta<nmaxBin; idelta++){
 276       //
 277       imin   =  TMath::Max(i-idelta,1);
 278       imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
 279       nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
 280       if (nsum==0) break;
 281       if (nsum>minEntries) break;
 282     }
 283     if (nsum<minEntries) continue;
 284     //
 285     hist->GetXaxis()->SetRange(imin,imax);
 286     projectionHist = hist->ProjectionY("gain",imin,imax);
 287     // Determine Median:
 288     Float_t xMin = projectionHist->GetXaxis()->GetXmin();
 289     Float_t xMax = projectionHist->GetXaxis()->GetXmax();
 290     Float_t xMedian = (xMin+xMax)*0.5;
 291     Float_t integral = 0;
 292     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
 293       integral+=projectionHist->GetBinContent(jbin);
 294     }
 295     //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
 296     //
 297     //
 298     Float_t currentSum=0;
 299     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
 300       currentSum += projectionHist->GetBinContent(jbin);
 301       if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
 302       if (currentSum<fracUp*integral)  xMax = projectionHist->GetBinCenter(jbin+1);
 303       if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
 304         xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
 305                    projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
 306           (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
 307       }
 308     }
 309     //
 310     Float_t rms  = projectionHist->GetRMS();
 311     //    i += interval;
 312     //
 313     Double_t xcenter =  hist->GetMean();
 314     Double_t xrms    =  hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
 315     Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
 316     if (rms>0){
 317       // cut on +- 4 RMS
 318       projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
 319       Double_t chi2 = funcGaus.GetChisquare();
 320       //
 321       xvec[counter] = xcenter;
 322       yvec[counter] = funcGaus.GetParameter(ival);
 323       xerr[counter] = xrms;
 324       yerr[counter] = funcGaus.GetParError(ival);
 325       if (useMedian) yvec[counter] = xMedian;
 326       if (cstream){
 327         (*cstream)<<"fitDebug"<<
 328           "xcenter="<<xcenter<<
 329           "xMin="<<xMin<<
 330           "xMax="<<xMax<<
 331           "xMedian="<<xMedian<<
 332           "xFitM"<<yvec[counter]<<
 333           "xFitS"<<yerr[counter]<<
 334           "chi2="<<chi2<<
 335           "\n";
 336       }
 337       counter++;
 338     }else{
 339       xvec[counter] = xcenter;
 340       yvec[counter] = xMedian;
 341       xerr[counter] = xrms;
 342       yerr[counter] = binWidth/TMath::Sqrt(12.);
 343       counter++;
 344     }
 345     delete projectionHist;
 346   }
 347
 348   TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
 349   delete [] xvec;
 350   delete [] yvec;
 351   delete [] xerr;
 352   delete [] yerr;
 353   delete hist;
 354   return graphErrors;
 355 }
 356
 357
 358 void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
 359
 360   // Method for the correct logarithmic binning of histograms
 361
 362   TAxis *axis = h->GetAxis(axisDim);
 363   int bins = axis->GetNbins();
 364
 365   Double_t from = axis->GetXmin();
 366   Double_t to = axis->GetXmax();
 367   Double_t *new_bins = new Double_t[bins + 1];
 368
 369   new_bins[0] = from;
 370   Double_t factor = pow(to/from, 1./bins);
 371
 372   for (int i = 1; i <= bins; i++) {
 373    new_bins[i] = factor * new_bins[i-1];
 374   }
 375   axis->Set(bins, new_bins);
 376   delete [] new_bins;
 377
 378 }
 379 void AliTPCcalibBase::BinLogX(TH1 *h) {
 380
 381   // Method for the correct logarithmic binning of histograms
 382
 383   TAxis *axis = h->GetXaxis();
 384   int bins = axis->GetNbins();
 385
 386   Double_t from = axis->GetXmin();
 387   Double_t to = axis->GetXmax();
 388   Double_t *new_bins = new Double_t[bins + 1];
 389
 390   new_bins[0] = from;
 391   Double_t factor = pow(to/from, 1./bins);
 392
 393   for (int i = 1; i <= bins; i++) {
 394    new_bins[i] = factor * new_bins[i-1];
 395   }
 396   axis->Set(bins, new_bins);
 397   delete [] new_bins;
 398
 399 }