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