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