CORRFW/AliCFTrackIsPrimaryCuts.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 // The class AliCFTrackIsPrimaryCut is designed to select reconstructed tracks
  17 // with a small impact parameter and tracks which are (not) daughters of kink
  18 // decays and to provide corresponding QA histograms.
  19 // This class inherits from the Analysis' Framework abstract base class
  20 // AliAnalysisCuts and is a part of the Correction Framework.
  21 // This class acts on single, reconstructed tracks, it is applicable on
  22 // ESD and AOD data.
  23 // It mainly consists of a IsSelected function that returns a boolean.
  24 // This function checks whether the considered track passes a set of cuts:
  25 // - min. and max. distance to main vertex in transverse plane (xy)
  26 // - min. and max. longitudinal distance to main vertex (z)
  27 // - min. and max. distance to main vertex as ellpise in xy - z plane
  28 // - all above cuts on absolute values or in units of sigma (resolution)
  29 // - min. and max. distance to main vertex in units of sigma (resolution)
  30 // - max. transverse (xy) and longitudinal (z) impact parameter resolution
  31 // - require that the dca calculation doesn't fail
  32 // - accept or not accept daughter tracks of kink decays
  33 //
  34 // By default, the distance to 'vertex calculated from tracks' is used.
  35 // Optionally the SPD (tracklet based) or TPC (TPC only tracks based) vertex
  36 // can be used.
  37 // Note: the distance to the TPC-vertex is already stored in the ESD,
  38 // the distance to the SPD-vertex has to be re-calculated by propagating each
  39 // track while executing this cut.
  40 //
  41 // The cut values for these cuts are set with the corresponding set functions.
  42 // All cut classes provided by the correction framework are supposed to be
  43 // added in the Analysis Framwork's class AliAnalysisFilter and applied by
  44 // the filter via a loop.
  45 //
  46 // author: I. Kraus (Ingrid.Kraus@cern.ch)
  47 // idea taken form
  48 // AliESDtrackCuts writte by Jan Fiete Grosse-Oetringhaus and
  49 // AliRsnDaughterCut class written by A. Pulvirenti.
  50
  51 #include <TCanvas.h>
  52 #include <TDirectory.h>
  53 #include <TH2.h>
  54 #include <TBits.h>
  55
  56 #include <AliESDtrack.h>
  57 #include <AliAODTrack.h>
  58 #include <AliESDEvent.h>
  59 #include <AliAODEvent.h>
  60 #include <AliLog.h>
  61 #include "AliCFTrackIsPrimaryCuts.h"
  62
  63 ClassImp(AliCFTrackIsPrimaryCuts)
  64
  65 //__________________________________________________________________________________
  66 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts() :
  67   AliCFCutBase(),
  68   fEvt(0x0),
  69   fUseSPDvertex(0),
  70   fUseTPCvertex(0),
  71   fMinDCAToVertexXY(0),
  72   fMinDCAToVertexZ(0),
  73   fMaxDCAToVertexXY(0),
  74   fMaxDCAToVertexZ(0),
  75   fDCAToVertex2D(0),
  76   fAbsDCAToVertex(0),
  77   fNSigmaToVertexMin(0),
  78   fNSigmaToVertexMax(0),
  79   fSigmaDCAxy(0),
  80   fSigmaDCAz(0),
  81   fRequireSigmaToVertex(0),
  82   fAODType(AliAODTrack::kUndef),
  83   fAcceptKinkDaughters(0),
  84   fhCutStatistics(0),
  85   fhCutCorrelation(0),
  86   fBitmap(0x0),
  87   fhNBinsNSigma(0),
  88   fhNBinsRequireSigma(0),
  89   fhNBinsAcceptKink(0),
  90   fhNBinsDcaXY(0),
  91   fhNBinsDcaZ(0),
  92   fhNBinsDcaXYnorm(0),
  93   fhNBinsDcaZnorm(0),
  94   fhNBinsSigmaDcaXY(0),
  95   fhNBinsSigmaDcaZ(0),
  96   fhBinLimNSigma(0x0),
  97   fhBinLimRequireSigma(0x0),
  98   fhBinLimAcceptKink(0x0),
  99   fhBinLimDcaXY(0x0),
 100   fhBinLimDcaZ(0x0),
 101   fhBinLimDcaXYnorm(0x0),
 102   fhBinLimDcaZnorm(0x0),
 103   fhBinLimSigmaDcaXY(0x0),
 104   fhBinLimSigmaDcaZ(0x0)
 105 {
 106   //
 107   // Default constructor
 108   //
 109   Initialise();
 110 }
 111 //__________________________________________________________________________________
 112 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(Char_t* name, Char_t* title) :
 113   AliCFCutBase(name,title),
 114   fEvt(0x0),
 115   fUseSPDvertex(0),
 116   fUseTPCvertex(0),
 117   fMinDCAToVertexXY(0),
 118   fMinDCAToVertexZ(0),
 119   fMaxDCAToVertexXY(0),
 120   fMaxDCAToVertexZ(0),
 121   fDCAToVertex2D(0),
 122   fAbsDCAToVertex(0),
 123   fNSigmaToVertexMin(0),
 124   fNSigmaToVertexMax(0),
 125   fSigmaDCAxy(0),
 126   fSigmaDCAz(0),
 127   fRequireSigmaToVertex(0),
 128   fAODType(AliAODTrack::kUndef),
 129   fAcceptKinkDaughters(0),
 130   fhCutStatistics(0),
 131   fhCutCorrelation(0),
 132   fBitmap(0x0),
 133   fhNBinsNSigma(0),
 134   fhNBinsRequireSigma(0),
 135   fhNBinsAcceptKink(0),
 136   fhNBinsDcaXY(0),
 137   fhNBinsDcaZ(0),
 138   fhNBinsDcaXYnorm(0),
 139   fhNBinsDcaZnorm(0),
 140   fhNBinsSigmaDcaXY(0),
 141   fhNBinsSigmaDcaZ(0),
 142   fhBinLimNSigma(0x0),
 143   fhBinLimRequireSigma(0x0),
 144   fhBinLimAcceptKink(0x0),
 145   fhBinLimDcaXY(0x0),
 146   fhBinLimDcaZ(0x0),
 147   fhBinLimDcaXYnorm(0x0),
 148   fhBinLimDcaZnorm(0x0),
 149   fhBinLimSigmaDcaXY(0x0),
 150   fhBinLimSigmaDcaZ(0x0)
 151 {
 152   //
 153   // Constructor
 154   //
 155   Initialise();
 156 }
 157 //__________________________________________________________________________________
 158 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(const AliCFTrackIsPrimaryCuts& c) :
 159   AliCFCutBase(c),
 160   fEvt(c.fEvt),
 161   fUseSPDvertex(c.fUseSPDvertex),
 162   fUseTPCvertex(c.fUseTPCvertex),
 163   fMinDCAToVertexXY(c.fMinDCAToVertexXY),
 164   fMinDCAToVertexZ(c.fMinDCAToVertexZ),
 165   fMaxDCAToVertexXY(c.fMaxDCAToVertexXY),
 166   fMaxDCAToVertexZ(c.fMaxDCAToVertexZ),
 167   fDCAToVertex2D(c.fDCAToVertex2D),
 168   fAbsDCAToVertex(c.fAbsDCAToVertex),
 169   fNSigmaToVertexMin(c.fNSigmaToVertexMin),
 170   fNSigmaToVertexMax(c.fNSigmaToVertexMax),
 171   fSigmaDCAxy(c.fSigmaDCAxy),
 172   fSigmaDCAz(c.fSigmaDCAz),
 173   fRequireSigmaToVertex(c.fRequireSigmaToVertex),
 174   fAODType(c.fAODType),
 175   fAcceptKinkDaughters(c.fAcceptKinkDaughters),
 176   fhCutStatistics(c.fhCutStatistics),
 177   fhCutCorrelation(c.fhCutCorrelation),
 178   fBitmap(c.fBitmap),
 179   fhNBinsNSigma(c.fhNBinsNSigma),
 180   fhNBinsRequireSigma(c.fhNBinsRequireSigma),
 181   fhNBinsAcceptKink(c.fhNBinsAcceptKink),
 182   fhNBinsDcaXY(c.fhNBinsDcaXY),
 183   fhNBinsDcaZ(c.fhNBinsDcaZ),
 184   fhNBinsDcaXYnorm(c.fhNBinsDcaXYnorm),
 185   fhNBinsDcaZnorm(c.fhNBinsDcaZnorm),
 186   fhNBinsSigmaDcaXY(c.fhNBinsSigmaDcaXY),
 187   fhNBinsSigmaDcaZ(c.fhNBinsSigmaDcaZ),
 188   fhBinLimNSigma(c.fhBinLimNSigma),
 189   fhBinLimRequireSigma(c.fhBinLimRequireSigma),
 190   fhBinLimAcceptKink(c.fhBinLimAcceptKink),
 191   fhBinLimDcaXY(c.fhBinLimDcaXY),
 192   fhBinLimDcaZ(c.fhBinLimDcaZ),
 193   fhBinLimDcaXYnorm(c.fhBinLimDcaXYnorm),
 194   fhBinLimDcaZnorm(c.fhBinLimDcaZnorm),
 195   fhBinLimSigmaDcaXY(c.fhBinLimSigmaDcaXY),
 196   fhBinLimSigmaDcaZ(c.fhBinLimSigmaDcaZ)
 197 {
 198   //
 199   // copy constructor
 200   //
 201   ((AliCFTrackIsPrimaryCuts &) c).Copy(*this);
 202 }
 203 //__________________________________________________________________________________
 204 AliCFTrackIsPrimaryCuts& AliCFTrackIsPrimaryCuts::operator=(const AliCFTrackIsPrimaryCuts& c)
 205 {
 206   //
 207   // Assignment operator
 208   //
 209   if (this != &c) {
 210     AliCFCutBase::operator=(c) ;
 211     fEvt = c.fEvt;
 212     fUseSPDvertex = c.fUseSPDvertex;
 213     fUseTPCvertex = c.fUseTPCvertex;
 214     fMinDCAToVertexXY = c.fMinDCAToVertexXY;
 215     fMinDCAToVertexZ = c.fMinDCAToVertexZ;
 216     fMaxDCAToVertexXY = c.fMaxDCAToVertexXY;
 217     fMaxDCAToVertexZ = c.fMaxDCAToVertexZ;
 218     fDCAToVertex2D = c.fDCAToVertex2D;
 219     fAbsDCAToVertex = c.fAbsDCAToVertex;
 220     fNSigmaToVertexMin = c.fNSigmaToVertexMin ;
 221     fNSigmaToVertexMax = c.fNSigmaToVertexMax ;
 222     fSigmaDCAxy = c.fSigmaDCAxy ;
 223     fSigmaDCAz = c.fSigmaDCAz ;
 224     fRequireSigmaToVertex = c.fRequireSigmaToVertex ;
 225     fAODType = c.fAODType ;
 226     fAcceptKinkDaughters = c.fAcceptKinkDaughters ;
 227     fhCutStatistics = c.fhCutStatistics ;
 228     fhCutCorrelation = c.fhCutCorrelation ;
 229     fBitmap =  c.fBitmap;
 230     fhNBinsNSigma = c.fhNBinsNSigma;
 231     fhNBinsRequireSigma = c.fhNBinsRequireSigma;
 232     fhNBinsAcceptKink = c.fhNBinsAcceptKink;
 233     fhNBinsDcaXY = c.fhNBinsDcaXY;
 234     fhNBinsDcaZ = c.fhNBinsDcaZ;
 235     fhNBinsDcaXYnorm = c.fhNBinsDcaXYnorm;
 236     fhNBinsDcaZnorm = c.fhNBinsDcaZnorm;
 237     fhNBinsSigmaDcaXY = c.fhNBinsSigmaDcaXY;
 238     fhNBinsSigmaDcaZ = c.fhNBinsSigmaDcaZ;
 239     fhBinLimNSigma = c.fhBinLimNSigma;
 240     fhBinLimRequireSigma = c.fhBinLimRequireSigma;
 241     fhBinLimAcceptKink = c.fhBinLimAcceptKink;
 242     fhBinLimDcaXY = c.fhBinLimDcaXY;
 243     fhBinLimDcaZ = c.fhBinLimDcaZ;
 244     fhBinLimDcaXYnorm = c.fhBinLimDcaXYnorm;
 245     fhBinLimDcaZnorm = c.fhBinLimDcaZnorm;
 246     fhBinLimSigmaDcaXY = c.fhBinLimSigmaDcaXY;
 247     fhBinLimSigmaDcaZ = c.fhBinLimSigmaDcaZ;
 248
 249     for (Int_t j=0; j<6; j++) fDCA[j] = c.fDCA[j];
 250     for (Int_t j=0; j<c.kNStepQA; j++){
 251       if(c.fhDcaXYvsDcaZ[j]) fhDcaXYvsDcaZ[j] = (TH2F*)c.fhDcaXYvsDcaZ[j]->Clone();
 252       for (Int_t i=0; i<c.kNHist; i++){
 253         if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
 254       }
 255     }
 256     ((AliCFTrackIsPrimaryCuts &) c).Copy(*this);
 257  }
 258   return *this;
 259 }
 260 //__________________________________________________________________________________
 261 AliCFTrackIsPrimaryCuts::~AliCFTrackIsPrimaryCuts()
 262 {
 263   //
 264   // destructor
 265   //
 266   if (fhCutStatistics)                  delete fhCutStatistics;
 267   if (fhCutCorrelation)                 delete fhCutCorrelation;
 268
 269   for (Int_t j=0; j<kNStepQA; j++){
 270     if(fhDcaXYvsDcaZ[j])        delete fhDcaXYvsDcaZ[j];
 271     for (Int_t i=0; i<kNHist; i++)
 272       if(fhQA[i][j])            delete fhQA[i][j];
 273   }
 274   if(fEvt)                      delete fEvt;
 275   if(fBitmap)                   delete fBitmap;
 276   if(fhBinLimNSigma)            delete fhBinLimNSigma;
 277   if(fhBinLimRequireSigma)      delete fhBinLimRequireSigma;
 278   if(fhBinLimAcceptKink)        delete fhBinLimAcceptKink;
 279   if(fhBinLimDcaXY)             delete fhBinLimDcaXY;
 280   if(fhBinLimDcaZ)              delete fhBinLimDcaZ;
 281   if(fhBinLimDcaXYnorm)         delete fhBinLimDcaXYnorm;
 282   if(fhBinLimDcaZnorm)          delete fhBinLimDcaZnorm;
 283   if(fhBinLimSigmaDcaXY)        delete fhBinLimSigmaDcaXY;
 284   if(fhBinLimSigmaDcaZ)         delete fhBinLimSigmaDcaZ;
 285 }
 286 //__________________________________________________________________________________
 287 void AliCFTrackIsPrimaryCuts::Initialise()
 288 {
 289   //
 290   // sets everything to zero
 291   //
 292   fUseSPDvertex = 0;
 293   fUseTPCvertex = 0;
 294   fMinDCAToVertexXY = 0;
 295   fMinDCAToVertexZ = 0;
 296   fMaxDCAToVertexXY = 0;
 297   fMaxDCAToVertexZ = 0;
 298   fDCAToVertex2D = 0;
 299   fAbsDCAToVertex = 0;
 300   fNSigmaToVertexMin = 0;
 301   fNSigmaToVertexMax = 0;
 302   fSigmaDCAxy = 0;
 303   fSigmaDCAz = 0;
 304   fRequireSigmaToVertex = 0;
 305   fAcceptKinkDaughters = 0;
 306   fAODType = AliAODTrack::kUndef;
 307
 308   SetMinDCAToVertexXY();
 309   SetMinDCAToVertexZ();
 310   SetMaxDCAToVertexXY();
 311   SetMaxDCAToVertexZ();
 312   SetDCAToVertex2D();
 313   SetAbsDCAToVertex();
 314   SetMinNSigmaToVertex();
 315   SetMaxNSigmaToVertex();
 316   SetMaxSigmaDCAxy();
 317   SetMaxSigmaDCAz();
 318   SetRequireSigmaToVertex();
 319   SetAcceptKinkDaughters();
 320   SetAODType();
 321
 322   for (Int_t j=0; j<6; j++) fDCA[j] = 0.;
 323   for (Int_t j=0; j<kNStepQA; j++)  {
 324     fhDcaXYvsDcaZ[j] = 0x0;
 325     for (Int_t i=0; i<kNHist; i++)
 326       fhQA[i][j] = 0x0;
 327   }
 328   fhCutStatistics = 0;
 329   fhCutCorrelation = 0;
 330   fBitmap=new TBits(0);
 331
 332   //set default bining for QA histograms
 333   SetHistogramBins(kCutNSigmaToVertex,100,0.,10.);
 334   SetHistogramBins(kCutRequireSigmaToVertex,5,-0.75,1.75);
 335   SetHistogramBins(kCutAcceptKinkDaughters,5,-0.75,1.75);
 336   SetHistogramBins(kDcaXY,500,-10.,10.);
 337   SetHistogramBins(kDcaZ,500,-10.,10.);
 338   SetHistogramBins(kDcaXYnorm,500,-10.,10.);
 339   SetHistogramBins(kDcaZnorm,500,-10.,10.);
 340   SetHistogramBins(kSigmaDcaXY,500,-0.1,0.9);
 341   SetHistogramBins(kSigmaDcaZ,500,-0.1,0.9);
 342 }
 343 //__________________________________________________________________________________
 344 void AliCFTrackIsPrimaryCuts::Copy(TObject &c) const
 345 {
 346   //
 347   // Copy function
 348   //
 349   AliCFTrackIsPrimaryCuts& target = (AliCFTrackIsPrimaryCuts &) c;
 350
 351   target.Initialise();
 352
 353   if (fhCutStatistics)  target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
 354   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
 355
 356   for (Int_t j=0; j<6; j++) target.fDCA[j] = fDCA[j];
 357   for (Int_t j=0; j<kNStepQA; j++){
 358     if(fhDcaXYvsDcaZ[j]) target.fhDcaXYvsDcaZ[j] = (TH2F*)fhDcaXYvsDcaZ[j]->Clone();
 359     for (Int_t i=0; i<kNHist; i++)
 360       if(fhQA[i][j]) target.fhQA[i][j] = (TH1F*)fhQA[i][j]->Clone();
 361   }
 362   TNamed::Copy(c);
 363 }
 364 //__________________________________________________________________________________
 365 void AliCFTrackIsPrimaryCuts::SetRecEventInfo(const TObject* evt) {
 366   //
 367   // Sets pointer to event information (AliESDEvent or AliAODEvent)
 368   //
 369   if (!evt) {
 370     AliError("Pointer to AliVEvent !");
 371     return;
 372   }
 373   TString className(evt->ClassName());
 374   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
 375     AliError("argument must point to an AliESDEvent or AliAODEvent !");
 376     return ;
 377   }
 378   fEvt = (AliVEvent*) evt;
 379 }
 380 //__________________________________________________________________________________
 381 void AliCFTrackIsPrimaryCuts::UseSPDvertex(Bool_t b) {
 382   fUseSPDvertex = b;
 383   if(fUseTPCvertex && fUseSPDvertex) {
 384         fUseSPDvertex = kFALSE;
 385         AliError("SPD and TPC vertex chosen. TPC vertex is preferred.");
 386   }
 387 }
 388 //__________________________________________________________________________________
 389 void AliCFTrackIsPrimaryCuts::UseTPCvertex(Bool_t b) {
 390   fUseTPCvertex = b;
 391   if(fUseTPCvertex) fUseSPDvertex = kFALSE;
 392 }
 393 //__________________________________________________________________________________
 394 void AliCFTrackIsPrimaryCuts::GetDCA(AliESDtrack* esdTrack)
 395 {
 396   if (!esdTrack) return;
 397
 398   Float_t b[2] = {0.,0.};
 399   Float_t bCov[3] = {0.,0.,0.};
 400   if(!fUseSPDvertex && !fUseTPCvertex)  esdTrack->GetImpactParameters(b,bCov);
 401   if( fUseTPCvertex)                    esdTrack->GetImpactParametersTPC(b,bCov);
 402
 403   if( fUseSPDvertex) {
 404         if (!fEvt) return;
 405         AliESDEvent * evt = 0x0 ;
 406         evt = dynamic_cast<AliESDEvent*>(fEvt);
 407         const AliESDVertex *vtx = evt->GetVertex();
 408         const Double_t Bz = evt->GetMagneticField();
 409         AliExternalTrackParam *cParam = 0;
 410         Bool_t success = esdTrack->RelateToVertex(vtx, Bz, kVeryBig, cParam);
 411         if (success) esdTrack->GetImpactParameters(b,bCov);
 412   }
 413
 414   if (bCov[0]<=0 || bCov[2]<=0) {
 415       bCov[0]=0; bCov[2]=0;
 416   }
 417   fDCA[0] = b[0]; // impact parameter xy
 418   fDCA[1] = b[1]; // impact parameter z
 419   fDCA[2] = TMath::Sqrt(bCov[0]); // resolution xy
 420   fDCA[3] = TMath::Sqrt(bCov[2]); // resolution z
 421
 422   if (!fAbsDCAToVertex) {
 423         if (fDCA[2] > 0) fDCA[0] = fDCA[0]/fDCA[2]; // normalised impact parameter xy
 424         if (fDCA[3] > 0) fDCA[1] = fDCA[1]/fDCA[3]; // normalised impact parameter z
 425   }
 426
 427   // get n_sigma
 428   if(!fUseSPDvertex && !fUseTPCvertex)
 429         fDCA[5] = AliESDtrackCuts::GetSigmaToVertex(esdTrack);
 430
 431   if(fUseTPCvertex) {
 432         fDCA[5] = -1;
 433         if (fDCA[2]==0 || fDCA[3]==0)
 434           return;
 435         fDCA[5] = 1000.;
 436         Float_t d = TMath::Sqrt(TMath::Power(b[0]/fDCA[2],2) + TMath::Power(b[1]/fDCA[3],2));
 437         if (TMath::Exp(-d * d / 2) < 1e-15)
 438           return;
 439         fDCA[5] = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 440   }
 441   return;
 442 }
 443 //__________________________________________________________________________________
 444 void AliCFTrackIsPrimaryCuts::GetDCA(AliAODTrack* aodTrack)
 445 {
 446   if (!aodTrack) return;
 447
 448   Double_t p[3] = {0.,0.,0.};
 449   Double_t cov[21] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 450
 451   aodTrack->XYZAtDCA(p); // position at DCA
 452   aodTrack->GetCovarianceXYZPxPyPz(cov);
 453
 454
 455   fDCA[5] = -1; // n_sigma
 456   if (p[0]==-999. || p[1]==-999. ||  p[2]==-999.) {
 457     AliError("dca info not available !");
 458     fDCA[0] = -999.; // impact parameter xy
 459     fDCA[1] = -999.; // impact parameter z
 460     return;
 461   }
 462
 463   if (!fEvt) return;
 464   AliAODEvent * evt = 0x0;
 465   evt = dynamic_cast<AliAODEvent*>(fEvt);
 466   // primary vertex is the "best": tracks, SPD or TPC vertex
 467   AliAODVertex * primaryVertex = evt->GetVertex(0);
 468   // dca = track postion - primary vertex position
 469   p[0] = p[0] - primaryVertex->GetX();
 470   p[1] = p[1] - primaryVertex->GetY();
 471   p[2] = p[2] - primaryVertex->GetZ();
 472   // calculate dca in transverse plane
 473   Float_t b[2] = {0.,0.};
 474   b[0] = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
 475   b[1] = p[2];
 476   // resolution
 477   Double_t bCov[3] = {0.,0.,0.};
 478   // how to calculate the resoultion in the transverse plane ?
 479   bCov[0] = 0.; // to do: calculate cov in transverse plane
 480   bCov[2] = 0.; // from cov in x and y, need to know correlation
 481
 482   if (bCov[0]<=0 || bCov[2]<=0) {
 483       bCov[0]=0; bCov[2]=0;
 484   }
 485   fDCA[0] = b[0]; // impact parameter xy
 486   fDCA[1] = b[1]; // impact parameter z
 487   fDCA[2] = TMath::Sqrt(bCov[0]); // resolution xy
 488   fDCA[3] = TMath::Sqrt(bCov[2]); // resolution z
 489
 490   if (!fAbsDCAToVertex) {
 491         AliError("resolution of impact parameter not available, use absolute dca cut instead !");
 492         if (fDCA[2] > 0) fDCA[0] = fDCA[0]/fDCA[2]; // normalised impact parameter xy
 493         if (fDCA[3] > 0) fDCA[1] = fDCA[1]/fDCA[3]; // normalised impact parameter z
 494   }
 495
 496   // get n_sigma
 497   if (fDCA[2]==0 || fDCA[3]==0)
 498         return;
 499   fDCA[5] = 1000.;
 500   Float_t d = TMath::Sqrt(TMath::Power(b[0]/fDCA[2],2) + TMath::Power(b[1]/fDCA[3],2));
 501   if (TMath::Exp(-d * d / 2) < 1e-15)
 502         return;
 503   fDCA[5] = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 504
 505   return;
 506 }
 507 //__________________________________________________________________________________
 508 void AliCFTrackIsPrimaryCuts::SelectionBitMap(TObject* obj)
 509 {
 510   //
 511   // test if the track passes the single cuts
 512   // and store the information in a bitmap
 513   //
 514
 515   // bitmap stores the decision of each single cut
 516   for(Int_t i=0; i<kNCuts; i++)fBitmap->SetBitNumber(i,kFALSE);
 517
 518   // check TObject and cast into ESDtrack
 519   if (!obj) return;
 520   if (!obj->InheritsFrom("AliVParticle")) {
 521     AliError("object must derived from AliVParticle !");
 522     return;
 523   }
 524
 525   Bool_t isESDTrack = strcmp(obj->ClassName(),"AliESDtrack") == 0 ? kTRUE : kFALSE ;
 526   Bool_t isAODTrack = strcmp(obj->ClassName(),"AliAODTrack") == 0 ? kTRUE : kFALSE ;
 527
 528   AliESDtrack * esdTrack = 0x0 ;
 529   AliAODTrack * aodTrack = 0x0 ;
 530   if (isESDTrack) esdTrack = dynamic_cast<AliESDtrack*>(obj);
 531   if (isAODTrack) aodTrack = dynamic_cast<AliAODTrack*>(obj);
 532
 533   // get the track to vertex parameter for ESD track
 534   if (isESDTrack) GetDCA(esdTrack);
 535   if (isAODTrack) GetDCA(aodTrack);
 536
 537   // check whether dca info is filled
 538   Bool_t dcaInfo = 0;
 539   if (fDCA[0]>-990. && fDCA[1]>-990.) dcaInfo = 1;
 540
 541   Float_t bxy = 0, bz = 0;
 542   bxy = TMath::Abs(fDCA[0]);
 543   bz  = TMath::Abs(fDCA[1]);
 544
 545   Float_t b2Dmin = 0, b2Dmax = 0;
 546   if (fMinDCAToVertexXY>0 && fMinDCAToVertexZ>0)
 547     b2Dmin = fDCA[0]*fDCA[0]/fMinDCAToVertexXY/fMinDCAToVertexXY + fDCA[1]*fDCA[1]/fMinDCAToVertexZ/fMinDCAToVertexZ;
 548   if (fMaxDCAToVertexXY>0 && fMaxDCAToVertexZ>0)
 549     b2Dmax = fDCA[0]*fDCA[0]/fMaxDCAToVertexXY/fMaxDCAToVertexXY + fDCA[1]*fDCA[1]/fMaxDCAToVertexZ/fMaxDCAToVertexZ;
 550
 551
 552   // fill the bitmap
 553   Int_t iCutBit = 0;
 554
 555   if (!dcaInfo || fDCAToVertex2D || (!fDCAToVertex2D && bxy >= fMinDCAToVertexXY && bxy <= fMaxDCAToVertexXY))
 556         fBitmap->SetBitNumber(iCutBit,kTRUE);
 557   iCutBit++;
 558
 559   if (!dcaInfo || fDCAToVertex2D || (!fDCAToVertex2D && bz  >= fMinDCAToVertexZ && bz  <= fMaxDCAToVertexZ))
 560         fBitmap->SetBitNumber(iCutBit,kTRUE);
 561   iCutBit++;
 562
 563   if (!dcaInfo || !fDCAToVertex2D || (fDCAToVertex2D && TMath::Sqrt(b2Dmin) > 1  && TMath::Sqrt(b2Dmax) < 1))
 564       fBitmap->SetBitNumber(iCutBit,kTRUE);
 565   iCutBit++;
 566
 567   if (!dcaInfo || (fDCA[5] >= fNSigmaToVertexMin && fDCA[5] <= fNSigmaToVertexMax))
 568       fBitmap->SetBitNumber(iCutBit,kTRUE);
 569   iCutBit++;
 570
 571   if (!dcaInfo || fDCA[2] < fSigmaDCAxy)
 572       fBitmap->SetBitNumber(iCutBit,kTRUE);
 573   iCutBit++;
 574
 575   if (!dcaInfo || fDCA[3] < fSigmaDCAz)
 576       fBitmap->SetBitNumber(iCutBit,kTRUE);
 577   iCutBit++;
 578
 579   if (!dcaInfo || !fRequireSigmaToVertex || (fDCA[5]>=0 && fRequireSigmaToVertex))
 580       fBitmap->SetBitNumber(iCutBit,kTRUE);
 581   iCutBit++;
 582
 583   if (!dcaInfo || fAcceptKinkDaughters || (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)<=0))
 584       fBitmap->SetBitNumber(iCutBit,kTRUE);
 585   iCutBit++;
 586
 587   if (isAODTrack) {
 588     if (fAODType==AliAODTrack::kUndef || fAODType == aodTrack->GetType()) {
 589       fBitmap->SetBitNumber(iCutBit,kTRUE);
 590     }
 591   }
 592   else fBitmap->SetBitNumber(iCutBit,kTRUE);
 593
 594   return;
 595 }
 596 //__________________________________________________________________________________
 597 Bool_t AliCFTrackIsPrimaryCuts::IsSelected(TObject* obj) {
 598   //
 599   // loops over decisions of single cuts and returns if the track is accepted
 600   //
 601   SelectionBitMap(obj);
 602
 603   if (fIsQAOn) FillHistograms(obj,0);
 604   Bool_t isSelected = kTRUE;
 605
 606   for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++) {
 607     if(!fBitmap->TestBitNumber(icut)) {
 608         isSelected = kFALSE;
 609         break;
 610     }
 611   }
 612
 613   if (!isSelected) return kFALSE ;
 614   if (fIsQAOn) FillHistograms(obj,1);
 615   return kTRUE;
 616 }
 617 //__________________________________________________________________________________
 618 void AliCFTrackIsPrimaryCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t *bins)
 619 {
 620   //
 621   // variable bin size
 622   //
 623   if(!fIsQAOn) return;
 624
 625   switch(index){
 626   case kCutNSigmaToVertex:
 627     fhNBinsNSigma=nbins+1;
 628     fhBinLimNSigma=new Double_t[nbins+1];
 629     for(Int_t i=0;i<nbins+1;i++)fhBinLimNSigma[i]=bins[i];
 630     break;
 631
 632   case kCutRequireSigmaToVertex:
 633     fhNBinsRequireSigma=nbins+1;
 634     fhBinLimRequireSigma=new Double_t[nbins+1];
 635     for(Int_t i=0;i<nbins+1;i++)fhBinLimRequireSigma[i]=bins[i];
 636     break;
 637
 638   case kCutAcceptKinkDaughters:
 639     fhNBinsAcceptKink=nbins+1;
 640     fhBinLimAcceptKink=new Double_t[nbins+1];
 641     for(Int_t i=0;i<nbins+1;i++)fhBinLimAcceptKink[i]=bins[i];
 642     break;
 643
 644   case kDcaXY:
 645     fhNBinsDcaXY=nbins+1;
 646     fhBinLimDcaXY=new Double_t[nbins+1];
 647     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXY[i]=bins[i];
 648     break;
 649
 650   case kDcaZ:
 651     fhNBinsDcaZ=nbins+1;
 652     fhBinLimDcaZ=new Double_t[nbins+1];
 653     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZ[i]=bins[i];
 654     break;
 655
 656   case kDcaXYnorm:
 657     fhNBinsDcaXYnorm=nbins+1;
 658     fhBinLimDcaXYnorm=new Double_t[nbins+1];
 659     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXYnorm[i]=bins[i];
 660     break;
 661
 662   case kDcaZnorm:
 663     fhNBinsDcaZnorm=nbins+1;
 664     fhBinLimDcaZnorm=new Double_t[nbins+1];
 665     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZnorm[i]=bins[i];
 666     break;
 667
 668   case kSigmaDcaXY:
 669     fhNBinsSigmaDcaXY=nbins+1;
 670     fhBinLimSigmaDcaXY=new Double_t[nbins+1];
 671     for(Int_t i=0;i<nbins+1;i++)fhBinLimSigmaDcaXY[i]=bins[i];
 672     break;
 673
 674   case kSigmaDcaZ:
 675     fhNBinsSigmaDcaZ=nbins+1;
 676     fhBinLimSigmaDcaZ=new Double_t[nbins+1];
 677     for(Int_t i=0;i<nbins+1;i++)fhBinLimSigmaDcaZ[i]=bins[i];
 678     break;
 679   }
 680 }
 681 //__________________________________________________________________________________
 682 void AliCFTrackIsPrimaryCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax)
 683 {
 684   //
 685   // fixed bin size
 686   //
 687   switch(index){
 688   case kCutNSigmaToVertex:
 689     fhNBinsNSigma=nbins+1;
 690     fhBinLimNSigma=new Double_t[nbins+1];
 691     for(Int_t i=0;i<nbins+1;i++)fhBinLimNSigma[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 692     break;
 693
 694   case kCutRequireSigmaToVertex:
 695     fhNBinsRequireSigma=nbins+1;
 696     fhBinLimRequireSigma=new Double_t[nbins+1];
 697     for(Int_t i=0;i<nbins+1;i++)fhBinLimRequireSigma[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 698     break;
 699
 700   case kCutAcceptKinkDaughters:
 701     fhNBinsAcceptKink=nbins+1;
 702     fhBinLimAcceptKink=new Double_t[nbins+1];
 703     for(Int_t i=0;i<nbins+1;i++)fhBinLimAcceptKink[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 704     break;
 705
 706   case kDcaXY:
 707     fhNBinsDcaXY=nbins+1;
 708     fhBinLimDcaXY=new Double_t[nbins+1];
 709     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXY[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 710     break;
 711
 712   case kDcaZ:
 713     fhNBinsDcaZ=nbins+1;
 714     fhBinLimDcaZ=new Double_t[nbins+1];
 715     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZ[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 716     break;
 717
 718   case kDcaXYnorm:
 719     fhNBinsDcaXYnorm=nbins+1;
 720     fhBinLimDcaXYnorm=new Double_t[nbins+1];
 721     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXYnorm[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 722     break;
 723
 724   case kDcaZnorm:
 725     fhNBinsDcaZnorm=nbins+1;
 726     fhBinLimDcaZnorm=new Double_t[nbins+1];
 727     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZnorm[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 728     break;
 729
 730   case kSigmaDcaXY:
 731     fhNBinsSigmaDcaXY=nbins+1;
 732     fhBinLimSigmaDcaXY=new Double_t[nbins+1];
 733     for(Int_t i=0;i<nbins+1;i++)fhBinLimSigmaDcaXY[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 734     break;
 735
 736   case kSigmaDcaZ:
 737     fhNBinsSigmaDcaZ=nbins+1;
 738     fhBinLimSigmaDcaZ=new Double_t[nbins+1];
 739     for(Int_t i=0;i<nbins+1;i++)fhBinLimSigmaDcaZ[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
 740     break;
 741   }
 742 }
 743 //__________________________________________________________________________________
 744  void AliCFTrackIsPrimaryCuts::DefineHistograms() {
 745   //
 746   // histograms for cut variables, cut statistics and cut correlations
 747   //
 748   Int_t color = 2;
 749
 750   // book cut statistics and cut correlation histograms
 751   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()), Form("%s cut statistics",GetName()), kNCuts,0.5,kNCuts+0.5);
 752   fhCutStatistics->SetLineWidth(2);
 753   fhCutStatistics->GetXaxis()->SetBinLabel(1,"dca xy");
 754   fhCutStatistics->GetXaxis()->SetBinLabel(2,"dca z");
 755   fhCutStatistics->GetXaxis()->SetBinLabel(3,"dca ellipse");
 756   fhCutStatistics->GetXaxis()->SetBinLabel(4,"n dca");
 757   fhCutStatistics->GetXaxis()->SetBinLabel(5,"sigma dca xy");
 758   fhCutStatistics->GetXaxis()->SetBinLabel(6,"sigma dca z");
 759   fhCutStatistics->GetXaxis()->SetBinLabel(7,"require dca");
 760   fhCutStatistics->GetXaxis()->SetBinLabel(8,"kink daughter");
 761   fhCutStatistics->GetXaxis()->SetBinLabel(9,"AOD type");
 762
 763   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()), Form("%s cut  correlation",GetName()), kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
 764   fhCutCorrelation->SetLineWidth(2);
 765   fhCutCorrelation->GetXaxis()->SetBinLabel(1,"dca xy");
 766   fhCutCorrelation->GetXaxis()->SetBinLabel(2,"dca z");
 767   fhCutCorrelation->GetXaxis()->SetBinLabel(3,"dca ellipse");
 768   fhCutCorrelation->GetXaxis()->SetBinLabel(4,"n dca");
 769   fhCutCorrelation->GetXaxis()->SetBinLabel(5,"sigma dca xy");
 770   fhCutCorrelation->GetXaxis()->SetBinLabel(6,"sigma dca z");
 771   fhCutCorrelation->GetXaxis()->SetBinLabel(7,"require dca");
 772   fhCutCorrelation->GetXaxis()->SetBinLabel(8,"kink daughter");
 773   fhCutCorrelation->GetXaxis()->SetBinLabel(9,"AOD type");
 774
 775   fhCutCorrelation->GetYaxis()->SetBinLabel(1,"dca xy");
 776   fhCutCorrelation->GetYaxis()->SetBinLabel(2,"dca z");
 777   fhCutCorrelation->GetYaxis()->SetBinLabel(3,"dca ellipse");
 778   fhCutCorrelation->GetYaxis()->SetBinLabel(4,"n dca");
 779   fhCutCorrelation->GetYaxis()->SetBinLabel(5,"sigma dca xy");
 780   fhCutCorrelation->GetYaxis()->SetBinLabel(6,"sigma dca z");
 781   fhCutCorrelation->GetYaxis()->SetBinLabel(7,"require dca");
 782   fhCutCorrelation->GetYaxis()->SetBinLabel(8,"kink daughter");
 783   fhCutCorrelation->GetYaxis()->SetBinLabel(9,"AOD type");
 784
 785   // book QA histograms
 786   Char_t str[256];
 787   for (Int_t i=0; i<kNStepQA; i++) {
 788     if (i==0) sprintf(str," ");
 789     else sprintf(str,"_cut");
 790
 791     fhDcaXYvsDcaZ[i]            = new  TH2F(Form("%s_dcaXYvsDcaZ%s",GetName(),str),"",200,-10,10,200,-10,10);
 792     fhQA[kCutNSigmaToVertex][i] = new TH1F(Form("%s_nSigmaToVertex%s",GetName(),str),"",fhNBinsNSigma-1,fhBinLimNSigma);
 793     fhQA[kCutRequireSigmaToVertex][i] = new TH1F(Form("%s_requireSigmaToVertex%s",GetName(),str),"",fhNBinsRequireSigma-1,fhBinLimRequireSigma);
 794     fhQA[kCutAcceptKinkDaughters][i] = new TH1F(Form("%s_acceptKinkDaughters%s",GetName(),str),"",fhNBinsAcceptKink-1,fhBinLimAcceptKink);
 795     fhQA[kDcaXY][i]             = new TH1F(Form("%s_dcaXY%s",GetName(),str),"",fhNBinsDcaXY-1,fhBinLimDcaXY);
 796     fhQA[kDcaZ][i]              = new TH1F(Form("%s_dcaZ%s",GetName(),str),"",fhNBinsDcaZ-1,fhBinLimDcaZ);
 797     fhQA[kDcaXYnorm][i]         = new TH1F(Form("%s_dcaXYnorm%s",GetName(),str),"",fhNBinsDcaXYnorm-1,fhBinLimDcaXYnorm);
 798     fhQA[kDcaZnorm][i]          = new TH1F(Form("%s_dcaZnorm%s",GetName(),str),"",fhNBinsDcaZnorm-1,fhBinLimDcaZnorm);
 799     fhQA[kSigmaDcaXY][i]        = new TH1F(Form("%s_sigmaDcaXY%s",GetName(),str),"",fhNBinsSigmaDcaXY-1,fhBinLimSigmaDcaXY);
 800     fhQA[kSigmaDcaZ][i]         = new TH1F(Form("%s_sigmaDcaZ%s",GetName(),str),"",fhNBinsSigmaDcaZ-1,fhBinLimSigmaDcaZ);
 801
 802     fhDcaXYvsDcaZ[i]->SetXTitle("impact par. d_{z}");
 803     fhDcaXYvsDcaZ[i]->SetYTitle("impact par. d_{xy}");
 804
 805     fhQA[kCutNSigmaToVertex][i]->SetXTitle("n #sigma to vertex");
 806     fhQA[kCutRequireSigmaToVertex][i]->SetXTitle("require #sigma to vertex");
 807     fhQA[kCutAcceptKinkDaughters][i]->SetXTitle("accept kink daughters");
 808     fhQA[kDcaXY][i]->SetXTitle("impact par. d_{xy}");
 809     fhQA[kDcaZ][i]->SetXTitle("impact par. d_{z}");
 810     fhQA[kDcaXYnorm][i]->SetXTitle("norm. impact par. d_{xy} / #sigma_{xy}");
 811     fhQA[kDcaZnorm][i]->SetXTitle("norm. impact par. d_{z} / #sigma_{z}");
 812     fhQA[kSigmaDcaXY][i]->SetXTitle("impact par. resolution #sigma_{xy}");
 813     fhQA[kSigmaDcaZ][i]->SetXTitle("impact par. resolution #sigma_{z}");
 814   }
 815   for(Int_t i=0; i<kNHist; i++) fhQA[i][1]->SetLineColor(color);
 816 }
 817 //__________________________________________________________________________________
 818 void AliCFTrackIsPrimaryCuts::FillHistograms(TObject* obj, Bool_t f)
 819 {
 820   //
 821   // fill the QA histograms
 822   //
 823
 824   if (!obj) return;
 825
 826   Bool_t isESDTrack = strcmp(obj->ClassName(),"AliESDtrack") == 0 ? kTRUE : kFALSE ;
 827   Bool_t isAODTrack = strcmp(obj->ClassName(),"AliAODTrack") == 0 ? kTRUE : kFALSE ;
 828
 829   AliESDtrack * esdTrack = 0x0 ;
 830   AliAODTrack * aodTrack = 0x0 ;
 831   if (isESDTrack) esdTrack = dynamic_cast<AliESDtrack*>(obj);
 832   if (isAODTrack) aodTrack = dynamic_cast<AliAODTrack*>(obj);
 833
 834   // f = 0: fill histograms before cuts
 835   // f = 1: fill histograms after cuts
 836
 837   // get the track to vertex parameter for ESD track
 838   if (esdTrack) {
 839
 840     fhQA[kDcaZ][f]->Fill(fDCA[1]);
 841     fhQA[kDcaXY][f]->Fill(fDCA[0]);
 842     fhDcaXYvsDcaZ[f]->Fill(fDCA[1],fDCA[0]);
 843     fhQA[kSigmaDcaXY][f]->Fill(fDCA[2]);
 844     fhQA[kSigmaDcaZ][f]->Fill(fDCA[3]);
 845 // // // // // // //  delete histograms
 846       fhQA[kDcaZnorm][f]->Fill(fDCA[1]);
 847       fhQA[kDcaXYnorm][f]->Fill(fDCA[0]);
 848
 849     fhQA[kCutNSigmaToVertex][f]->Fill(fDCA[5]);
 850     if (fDCA[5]<0 && fRequireSigmaToVertex)
 851       fhQA[kCutRequireSigmaToVertex][f]->Fill(0.);
 852     else
 853       fhQA[kCutRequireSigmaToVertex][f]->Fill(1.);
 854
 855     if (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
 856       fhQA[kCutAcceptKinkDaughters][f]->Fill(0.);
 857     else
 858       fhQA[kCutAcceptKinkDaughters][f]->Fill(1.);
 859   }
 860
 861   // fill cut statistics and cut correlation histograms with information from the bitmap
 862   if (f) return;
 863
 864   // Get the bitmap of the single cuts
 865   if ( !obj ) return;
 866   SelectionBitMap(obj);
 867
 868   // Number of single cuts in this class
 869   UInt_t ncuts = fBitmap->GetNbits();
 870   for(UInt_t bit=0; bit<ncuts;bit++) {
 871     if (!fBitmap->TestBitNumber(bit)) {
 872         fhCutStatistics->Fill(bit+1);
 873         for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
 874           if (!fBitmap->TestBitNumber(bit2))
 875             fhCutCorrelation->Fill(bit+1,bit2+1);
 876         }
 877     }
 878   }
 879 }
 880 //__________________________________________________________________________________
 881 void AliCFTrackIsPrimaryCuts::SaveHistograms(const Char_t* dir) {
 882   //
 883   // saves the histograms in a directory (dir)
 884   //
 885   if(!fIsQAOn) return;
 886
 887   if (!dir)
 888     dir = GetName();
 889
 890   gDirectory->mkdir(dir);
 891   gDirectory->cd(dir);
 892
 893   gDirectory->mkdir("before_cuts");
 894   gDirectory->mkdir("after_cuts");
 895
 896   fhCutStatistics->Write();
 897   fhCutCorrelation->Write();
 898
 899   for (Int_t j=0; j<kNStepQA; j++) {
 900     if (j==0)
 901       gDirectory->cd("before_cuts");
 902     else
 903       gDirectory->cd("after_cuts");
 904
 905     fhDcaXYvsDcaZ[j]    ->Write();
 906
 907     for(Int_t i=0; i<kNHist; i++) fhQA[i][j]->Write();
 908
 909     gDirectory->cd("../");
 910   }
 911   gDirectory->cd("../");
 912 }
 913 //__________________________________________________________________________________
 914 void AliCFTrackIsPrimaryCuts::DrawHistograms()
 915 {
 916   //
 917   // draws some histograms
 918   //
 919   if(!fIsQAOn) return;
 920
 921   // pad margins
 922   Float_t right = 0.03;
 923   Float_t left = 0.175;
 924   Float_t top = 0.03;
 925   Float_t bottom = 0.175;
 926
 927   TCanvas* canvas1 = new TCanvas("Track_QA_Primary_1", "Track QA Primary 1", 800, 500);
 928   canvas1->Divide(2, 1);
 929
 930   canvas1->cd(1);
 931   fhCutStatistics->SetStats(kFALSE);
 932   fhCutStatistics->LabelsOption("v");
 933   gPad->SetLeftMargin(left);
 934   gPad->SetBottomMargin(0.25);
 935   gPad->SetRightMargin(right);
 936   gPad->SetTopMargin(0.1);
 937   fhCutStatistics->Draw();
 938
 939   canvas1->cd(2);
 940   fhCutCorrelation->SetStats(kFALSE);
 941   fhCutCorrelation->LabelsOption("v");
 942   gPad->SetLeftMargin(0.30);
 943   gPad->SetRightMargin(bottom);
 944   gPad->SetTopMargin(0.1);
 945   gPad->SetBottomMargin(0.25);
 946   fhCutCorrelation->Draw("COLZ");
 947
 948   canvas1->SaveAs(Form("%s.eps", canvas1->GetName()));
 949   canvas1->SaveAs(Form("%s.ps", canvas1->GetName()));
 950
 951   // -----
 952
 953   TCanvas* canvas2 = new TCanvas("Track_QA_Primary_2", "Track QA Primary 2", 800, 800);
 954   canvas2->Divide(2, 2);
 955
 956   canvas2->cd(1);
 957   gPad->SetRightMargin(right);
 958   gPad->SetLeftMargin(left);
 959   gPad->SetTopMargin(top);
 960   gPad->SetBottomMargin(bottom);
 961   fhQA[kDcaXY][0]->SetStats(kFALSE);
 962   fhQA[kDcaXY][0]->Draw();
 963   fhQA[kDcaXY][1]->Draw("same");
 964
 965   canvas2->cd(2);
 966   gPad->SetRightMargin(right);
 967   gPad->SetLeftMargin(left);
 968   gPad->SetTopMargin(top);
 969   gPad->SetBottomMargin(bottom);
 970   fhQA[kDcaZ][0]->SetStats(kFALSE);
 971   fhQA[kDcaZ][0]->Draw();
 972   fhQA[kDcaZ][1]->Draw("same");
 973
 974   canvas2->cd(3);
 975 //   fhDXYvsDZ[0]->SetStats(kFALSE);
 976   gPad->SetLogz();
 977   gPad->SetLeftMargin(bottom);
 978   gPad->SetTopMargin(0.1);
 979   gPad->SetBottomMargin(bottom);
 980   gPad->SetRightMargin(0.2);
 981   fhDcaXYvsDcaZ[0]->Draw("COLZ");
 982
 983   canvas2->cd(4);
 984 //   fhDXYvsDZ[1]->SetStats(kFALSE);
 985   gPad->SetLogz();
 986   gPad->SetLeftMargin(bottom);
 987   gPad->SetTopMargin(0.1);
 988   gPad->SetBottomMargin(bottom);
 989   gPad->SetRightMargin(0.2);
 990   fhDcaXYvsDcaZ[1]->Draw("COLZ");
 991
 992   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
 993   canvas2->SaveAs(Form("%s.ps", canvas2->GetName()));
 994
 995   // -----
 996
 997   TCanvas* canvas3 = new TCanvas("Track_QA_Primary_3", "Track QA Primary 3", 800, 400);
 998   canvas3->Divide(2, 1);
 999
1000   canvas3->cd(1);
1001   gPad->SetRightMargin(right);
1002   gPad->SetLeftMargin(left);
1003   gPad->SetTopMargin(top);
1004   gPad->SetBottomMargin(bottom);
1005   fhQA[kDcaXYnorm][0]->SetStats(kFALSE);
1006   fhQA[kDcaXYnorm][0]->Draw();
1007   fhQA[kDcaXYnorm][1]->Draw("same");
1008
1009   canvas3->cd(2);
1010   gPad->SetRightMargin(right);
1011   gPad->SetLeftMargin(left);
1012   gPad->SetTopMargin(top);
1013   gPad->SetBottomMargin(bottom);
1014   fhQA[kDcaZnorm][0]->SetStats(kFALSE);
1015   fhQA[kDcaZnorm][0]->Draw();
1016   fhQA[kDcaZnorm][1]->Draw("same");
1017
1018   canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
1019   canvas3->SaveAs(Form("%s.ps", canvas3->GetName()));
1020
1021   // -----
1022
1023   TCanvas* canvas4 = new TCanvas("Track_QA_Primary_4", "Track QA Primary 4", 1200, 800);
1024   canvas4->Divide(3, 2);
1025
1026   canvas4->cd(1);
1027   gPad->SetRightMargin(right);
1028   gPad->SetLeftMargin(left);
1029   gPad->SetTopMargin(top);
1030   gPad->SetBottomMargin(bottom);
1031   fhQA[kSigmaDcaXY][0]->SetStats(kFALSE);
1032   fhQA[kSigmaDcaXY][0]->Draw();
1033   fhQA[kSigmaDcaXY][1]->Draw("same");
1034
1035   canvas4->cd(2);
1036   gPad->SetRightMargin(right);
1037   gPad->SetLeftMargin(left);
1038   gPad->SetTopMargin(top);
1039   gPad->SetBottomMargin(bottom);
1040   fhQA[kSigmaDcaZ][0]->SetStats(kFALSE);
1041   fhQA[kSigmaDcaZ][0]->Draw();
1042   fhQA[kSigmaDcaZ][1]->Draw("same");
1043
1044   canvas4->cd(4);
1045   gPad->SetRightMargin(right);
1046   gPad->SetLeftMargin(left);
1047   gPad->SetTopMargin(top);
1048   gPad->SetBottomMargin(bottom);
1049   fhQA[kCutNSigmaToVertex][0]->SetStats(kFALSE);
1050   fhQA[kCutNSigmaToVertex][0]->Draw();
1051   fhQA[kCutNSigmaToVertex][1]->Draw("same");
1052
1053   canvas4->cd(5);
1054   gPad->SetRightMargin(right);
1055   gPad->SetLeftMargin(left);
1056   gPad->SetTopMargin(top);
1057   gPad->SetBottomMargin(bottom);
1058   fhQA[kCutRequireSigmaToVertex][0]->SetStats(kFALSE);
1059   fhQA[kCutRequireSigmaToVertex][0]->Draw();
1060   fhQA[kCutRequireSigmaToVertex][1]->Draw("same");
1061
1062   canvas4->cd(6);
1063   gPad->SetRightMargin(right);
1064   gPad->SetLeftMargin(left);
1065   gPad->SetTopMargin(top);
1066   gPad->SetBottomMargin(bottom);
1067   fhQA[kCutAcceptKinkDaughters][0]->SetStats(kFALSE);
1068   fhQA[kCutAcceptKinkDaughters][0]->Draw();
1069   fhQA[kCutAcceptKinkDaughters][1]->Draw("same");
1070
1071   canvas4->SaveAs(Form("%s.eps", canvas4->GetName()));
1072   canvas4->SaveAs(Form("%s.ps", canvas4->GetName()));
1073 }
1074 //__________________________________________________________________________________
1075 void AliCFTrackIsPrimaryCuts::AddQAHistograms(TList *qaList) {
1076   //
1077   // saves the histograms in a TList
1078   //
1079   DefineHistograms();
1080
1081   qaList->Add(fhCutStatistics);
1082   qaList->Add(fhCutCorrelation);
1083
1084   for (Int_t j=0; j<kNStepQA; j++) {
1085     qaList->Add(fhDcaXYvsDcaZ[j]);
1086     for(Int_t i=0; i<kNHist; i++)
1087         qaList->Add(fhQA[i][j]);
1088   }
1089 }