PWGHF/hfe/AliHFEextraCuts.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 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
  17 // Cuts stored here:
  18 // - ITS pixels
  19 // - TPC cluster ratio
  20 // - TRD tracklets
  21 //
  22 // Authors:
  23 //   Markus Fasel <M.Fasel@gsi.de>
  24 //
  25 #include <TBits.h>
  26 #include <TClass.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TList.h>
  30 #include <TString.h>
  31 #include <TMath.h>
  32
  33 #include "AliAODTrack.h"
  34 #include "AliAODPid.h"
  35 #include "AliAODVertex.h"
  36 #include "AliAODTrack.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDVertex.h"
  39 #include "AliESDtrack.h"
  40 #include "AliLog.h"
  41 #include "AliMCParticle.h"
  42 #include "AliVEvent.h"
  43 #include "AliVTrack.h"
  44 #include "AliVParticle.h"
  45 #include "AliVertexerTracks.h"
  46 #include "AliVVertex.h"
  47
  48 #include "AliHFEextraCuts.h"
  49
  50 ClassImp(AliHFEextraCuts)
  51
  52 const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
  53
  54 //______________________________________________________
  55 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
  56   AliCFCutBase(name, title),
  57   fEvent(NULL),
  58   fCutCorrelation(0),
  59   fRequirements(0),
  60   fMinNClustersTPC(0),
  61   fMinNClustersTPCPID(0),
  62   fClusterRatioTPC(0.),
  63   fMinTrackletsTRD(0),
  64   fMinNbITScls(0),
  65   fTRDtrackletsExact(0),
  66   fPixelITS(0),
  67   fTPCclusterDef(0),
  68   fTPCclusterRatioDef(0),
  69   fFractionTPCShared(-1.0),
  70   fAbsHFEImpactParamNsigmaR(kTRUE),
  71   fCheck(kFALSE),
  72   fQAlist(0x0) ,
  73   fDebugLevel(0)
  74 {
  75   //
  76   // Default Constructor
  77   //
  78   //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
  79   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
  80   memset(fIPcutParam, 0, sizeof(Float_t) * 4);
  81 }
  82
  83 //______________________________________________________
  84 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
  85   AliCFCutBase(c),
  86   fEvent(c.fEvent),
  87   fCutCorrelation(c.fCutCorrelation),
  88   fRequirements(c.fRequirements),
  89   fMinNClustersTPC(c.fMinNClustersTPC),
  90   fMinNClustersTPCPID(c.fMinNClustersTPCPID),
  91   fClusterRatioTPC(c.fClusterRatioTPC),
  92   fMinTrackletsTRD(c.fMinTrackletsTRD),
  93   fMinNbITScls(c.fMinNbITScls),
  94   fTRDtrackletsExact(c.fTRDtrackletsExact),
  95   fPixelITS(c.fPixelITS),
  96   fTPCclusterDef(c.fTPCclusterDef),
  97   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
  98   fFractionTPCShared(c.fFractionTPCShared),
  99   fAbsHFEImpactParamNsigmaR(c.fAbsHFEImpactParamNsigmaR),
 100   fCheck(c.fCheck),
 101   fQAlist(0x0),
 102   fDebugLevel(0)
 103 {
 104   //
 105   // Copy constructor
 106   // Performs a deep copy
 107   //
 108   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 109   memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 110   if(IsQAOn()){
 111     fIsQAOn = kTRUE;
 112     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 113     if(fQAlist) fQAlist->SetOwner();
 114   }
 115 }
 116
 117 //______________________________________________________
 118 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
 119   //
 120   // Assignment operator
 121   //
 122   if(this != &c){
 123     AliCFCutBase::operator=(c);
 124     fEvent = c.fEvent;
 125     fCutCorrelation = c.fCutCorrelation;
 126     fRequirements = c.fRequirements;
 127     fClusterRatioTPC = c.fClusterRatioTPC;
 128     fMinNClustersTPC = c.fMinNClustersTPC;
 129     fMinNClustersTPCPID = c.fMinNClustersTPCPID;
 130     fMinTrackletsTRD = c.fMinTrackletsTRD;
 131     fMinNbITScls = c.fMinNbITScls;
 132     fTRDtrackletsExact = c.fTRDtrackletsExact;
 133     fPixelITS = c.fPixelITS;
 134     fTPCclusterDef = c.fTPCclusterDef;
 135     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
 136     fFractionTPCShared = c.fFractionTPCShared;
 137     fAbsHFEImpactParamNsigmaR = c.fAbsHFEImpactParamNsigmaR;
 138     fCheck = c.fCheck;
 139     fDebugLevel = c.fDebugLevel;
 140     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 141     memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 142     if(IsQAOn()){
 143       fIsQAOn = kTRUE;
 144       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 145       if(fQAlist) fQAlist->SetOwner();
 146     }else fQAlist = 0x0;
 147   }
 148   return *this;
 149 }
 150
 151 //______________________________________________________
 152 AliHFEextraCuts::~AliHFEextraCuts(){
 153   //
 154   // Destructor
 155   //
 156   if(fQAlist) delete fQAlist;
 157 }
 158
 159 //______________________________________________________
 160 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
 161   //
 162   // Set Virtual event an make a copy
 163   //
 164   if (!event) {
 165     AliError("Pointer to AliVEvent !");
 166     return;
 167   }
 168   TString className(event->ClassName());
 169   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
 170     AliError("argument must point to an AliESDEvent or AliAODEvent !");
 171     return ;
 172   }
 173   fEvent = (AliVEvent*) event;
 174
 175 }
 176
 177 //______________________________________________________
 178 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
 179   //
 180   // Steering function for the track selection
 181   //
 182   TClass *type = o->IsA();
 183   AliDebug(2, Form("Object type %s", type->GetName()));
 184   if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
 185     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
 186   }
 187   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
 188 }
 189
 190 //______________________________________________________
 191 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
 192   //
 193   // Checks cuts on reconstructed tracks
 194   // returns true if track is selected
 195   // QA histograms are filled before track selection and for
 196   // selected tracks after track selection
 197   //
 198   AliDebug(1, "Called");
 199   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
 200   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
 201   // Apply cuts
 202   Float_t impactR, impactZ;
 203   Double_t hfeimpactR, hfeimpactnsigmaR;
 204   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
 205   Double_t maximpactRcut;
 206   GetImpactParameters(track, impactR, impactZ);
 207   if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 208     // Protection for PbPb
 209     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
 210     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
 211   }
 212   Int_t  nclsITS = GetITSNbOfcls(track);
 213   UInt_t nclsTPC = GetTPCncls(track);
 214   UInt_t nclsTPCPID = track->GetTPCsignalN();
 215   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 216   Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
 217   Double_t ratioTPC = GetTPCclusterRatio(track);
 218   UChar_t trdTracklets;
 219   trdTracklets = track->GetTRDntrackletsPID();
 220   UChar_t itsPixel = track->GetITSClusterMap();
 221   Int_t status1 = GetITSstatus(track, 0);
 222   Int_t status2 = GetITSstatus(track, 1);
 223   Bool_t statusL0 = CheckITSstatus(status1);
 224   Bool_t statusL1 = CheckITSstatus(status2);
 225   if(TESTBIT(fRequirements, kTPCfractionShared)) {
 226     // cut on max fraction of shared TPC clusters
 227     if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
 228   }
 229   if(TESTBIT(fRequirements, kMinImpactParamR)){
 230     // cut on min. Impact Parameter in Radial direction
 231     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
 232   }
 233   if(TESTBIT(fRequirements, kMinImpactParamZ)){
 234     // cut on min. Impact Parameter in Z direction
 235     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
 236   }
 237   if(TESTBIT(fRequirements, kMaxImpactParamR)){
 238     // cut on max. Impact Parameter in Radial direction
 239     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
 240   }
 241   if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
 242     // HFE Impact parameter cut
 243     GetMaxImpactParameterCutR(track,maximpactRcut);
 244     if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
 245   }
 246   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
 247     // cut on max. Impact Parameter in Z direction
 248     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
 249   }
 250   if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
 251     // cut on min. HFE Impact Parameter in Radial direction
 252     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
 253   }
 254   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 255     // cut on max. HFE Impact Parameter n sigma in Radial direction
 256     if(fAbsHFEImpactParamNsigmaR) {
 257 //      if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
 258       if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
 259         SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 260 //        printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
 261       }
 262     }
 263     else {
 264       if(hfeimpactnsigmaRcut>0){
 265         if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
 266           SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 267           //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
 268         }
 269       }
 270       else{
 271         if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
 272           SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 273           //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
 274         }
 275       }
 276     }
 277   }
 278   if(TESTBIT(fRequirements, kClusterRatioTPC)){
 279     // cut on min ratio of found TPC clusters vs findable TPC clusters
 280     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
 281   }
 282   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
 283     // cut on minimum number of TRD tracklets
 284     AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
 285     if(fTRDtrackletsExact){
 286       if(trdTracklets == fMinTrackletsTRD) {
 287         SETBIT(survivedCut, kMinTrackletsTRD);
 288         AliDebug(1, "Track Selected");
 289       }
 290     }else{
 291       if(trdTracklets >= fMinTrackletsTRD){
 292         SETBIT(survivedCut, kMinTrackletsTRD);
 293         AliDebug(1, "Track Selected");
 294       }
 295       //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
 296     }
 297   }
 298
 299   if(TESTBIT(fRequirements, kMinNbITScls)){
 300     // cut on minimum number of ITS clusters
 301     //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
 302     AliDebug(1, Form("Min ITS clusters: [%d|%d]\n", fMinNbITScls, nclsITS));
 303     if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
 304   }
 305
 306   if(TESTBIT(fRequirements, kMinNClustersTPC)){
 307     // cut on minimum number of TPC tracklets
 308     //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
 309     AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
 310     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
 311   }
 312   if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
 313     AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
 314     if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
 315   }
 316   if(TESTBIT(fRequirements, kPixelITS)){
 317     // cut on ITS pixel layers
 318     AliDebug(1, "ITS cluster Map: ");
 319     //PrintBitMap(itsPixel);
 320     switch(fPixelITS){
 321       case kFirst:
 322         AliDebug(2, "First");
 323         //printf("First\n");
 324               if(TESTBIT(itsPixel, 0))
 325                 SETBIT(survivedCut, kPixelITS);
 326         else
 327                 if(fCheck && !statusL0)
 328                   SETBIT(survivedCut, kPixelITS);
 329                     break;
 330       case kSecond:
 331         //printf("Second\n");
 332         AliDebug(2, "Second");
 333               if(TESTBIT(itsPixel, 1))
 334                 SETBIT(survivedCut, kPixelITS);
 335         else
 336                 if(fCheck && !statusL1)
 337                   SETBIT(survivedCut, kPixelITS);
 338                     break;
 339       case kBoth:
 340         //printf("Both\n");
 341         AliDebug(2, "Both");
 342               if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
 343                       SETBIT(survivedCut, kPixelITS);
 344               else
 345           if(fCheck && !(statusL0 && statusL1))
 346                         SETBIT(survivedCut, kPixelITS);
 347               break;
 348       case kAny:
 349         //printf("Any\n");
 350         AliDebug(2, "Any");
 351               if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
 352                 SETBIT(survivedCut, kPixelITS);
 353         else
 354                 if(fCheck && !(statusL0 || statusL1))
 355                     SETBIT(survivedCut, kPixelITS);
 356                     break;
 357       case kExclusiveSecond:
 358         AliDebug(2, "Exlusive second");
 359         if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
 360           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
 361             SETBIT(survivedCut, kPixelITS);
 362         } else {
 363           if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
 364             SETBIT(survivedCut, kPixelITS);
 365         }
 366         break;
 367       case kNone:
 368         // No cut applied, set as survived
 369         SETBIT(survivedCut, kPixelITS);
 370         break;
 371       default:
 372         // default, defined as no cut applied
 373         AliDebug(2, "None");
 374         SETBIT(survivedCut, kPixelITS);
 375         break;
 376     }
 377     AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
 378   }
 379
 380   if(TESTBIT(fRequirements, kTOFPID)){
 381     // cut on TOF pid
 382     if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
 383   }
 384
 385   if(TESTBIT(fRequirements, kTOFmismatch)){
 386     // cut on TOF mismatch
 387     if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
 388   }
 389
 390   if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
 391       // cut on TPC PID cleanup
 392       Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
 393       if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
 394   }
 395
 396   if(TESTBIT(fRequirements, kEMCALmatch)){
 397     if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
 398   }
 399
 400   if(TESTBIT(fRequirements, kRejectKinkDaughter)){
 401     //printf("test daughter\n");
 402     if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
 403     //else printf("Found kink daughter\n");
 404   }
 405
 406   if(TESTBIT(fRequirements, kRejectKinkMother)){
 407     //printf("test mother\n");
 408     if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
 409     //else printf("Found kink mother\n");
 410   }
 411
 412   if(fRequirements == survivedCut){
 413     //
 414     // Track selected
 415     //
 416     AliDebug(2, "Track Survived cuts\n");
 417     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
 418     return kTRUE;
 419   }
 420   AliDebug(2, "Track cut");
 421   if(IsQAOn()) FillCutCorrelation(survivedCut);
 422   return kFALSE;
 423 }
 424
 425 //______________________________________________________
 426 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
 427   //
 428   // Checks cuts on Monte Carlo tracks
 429   // returns true if track is selected
 430   // QA histograms are filled before track selection and for
 431   // selected tracks after track selection
 432   //
 433   return kTRUE; // not yet implemented
 434 }
 435
 436 //______________________________________________________
 437 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
 438   //
 439   // Fill the QA histograms for ESD tracks
 440   // Function can be called before cuts or after cut application (second argument)
 441   //
 442   Float_t impactR, impactZ;
 443   GetImpactParameters(track, impactR, impactZ);
 444   TH1 *htmp = NULL;
 445   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
 446   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
 447   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 448   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
 449   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
 450   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
 451   UChar_t itsPixel = track->GetITSClusterMap();
 452   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
 453   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 454   if(pixelHist){
 455     Double_t firstEntry = 0.5;
 456     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 457       pixelHist->Fill(firstEntry + 3);
 458     else{
 459       if(itsPixel & BIT(0)){
 460         pixelHist->Fill(firstEntry);
 461         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 462         else pixelHist->Fill(firstEntry + 4);
 463       }
 464       if(itsPixel & BIT(1)){
 465         pixelHist->Fill(firstEntry + 1);
 466         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 467       }
 468     }
 469   }
 470   // Fill histogram with the status bits
 471   TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
 472   if(hStatusBits) {
 473     hStatusBits->Fill(0);  // Fill first bin with all tracks
 474     if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
 475     if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
 476     if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
 477     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
 478   }
 479   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
 480 }
 481
 482 // //______________________________________________________
 483 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 484 //   //
 485 //   // Fill the QA histograms for MC tracks
 486 //   // Function can be called before cuts or after cut application (second argument)
 487 //   // Not yet implemented
 488 //   //
 489 // }
 490
 491 //______________________________________________________
 492 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 493   //
 494   // Fill cut correlation histograms for tracks that didn't pass cuts
 495   //
 496   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
 497   if(correlation){
 498     for(Int_t icut = 0; icut < kNcuts; icut++){
 499       if(!TESTBIT(fRequirements, icut)) continue;
 500       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 501         if(!TESTBIT(fRequirements, jcut)) continue;
 502         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 503                 correlation->Fill(icut, jcut);
 504       }
 505     }
 506   }
 507 }
 508
 509 //______________________________________________________
 510 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 511   //
 512   // Add QA histograms
 513   // For each cut a histogram before and after track cut is created
 514   // Histos before respectively after cut are stored in different lists
 515   // Additionally a histogram with the cut correlation is created and stored
 516   // in the top directory
 517   //
 518
 519   TH1 *histo1D = 0x0;
 520   TH2 *histo2D = 0x0;
 521   TString cutstr[2] = {"before", "after"};
 522
 523   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
 524   for(Int_t icond = 0; icond < 2; icond++){
 525     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
 526     fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
 527     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 528     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 529     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
 530     fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
 531     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 532     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 533     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
 534     fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
 535     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 536     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 537     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
 538     fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
 539     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 540     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 541     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
 542     fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
 543     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
 544     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 545     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
 546     fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
 547     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 548     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 549     Int_t first = histo1D->GetXaxis()->GetFirst();
 550     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 551     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 552       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 553     qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
 554     fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
 555     histo1D->GetXaxis()->SetTitle("Track Status Bit");
 556     histo1D->GetYaxis()->SetTitle("Number of tracks");
 557     TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
 558     for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
 559     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcdEdxClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters for dEdx calculation", 161, 0, 160)), 7 + icond * fgkNQAhistos);
 560     fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
 561     histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
 562     histo1D->GetYaxis()->SetTitle("counts");
 563   }
 564   // Add cut correlation
 565   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
 566   fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
 567   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters" "Fraction Shared TPC clusters", "TOFPID", "No TOFmismatch", "EMCALmatch", "ImpactParam"};
 568   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 569   for(Int_t icut = 0; icut < kNcuts; icut++){
 570     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 571     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 572   }
 573 }
 574
 575 //______________________________________________________
 576 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 577   for(Int_t ibit = 32; ibit--; )
 578     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 579   printf("\n");
 580 }
 581
 582 //______________________________________________________
 583 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 584   //
 585   // Check whether ITS area is dead
 586   //
 587   Bool_t status;
 588   switch(itsStatus){
 589     case 2: status = kFALSE; break;
 590     case 3: status = kFALSE; break;
 591     case 7: status = kFALSE; break;
 592     default: status = kTRUE;
 593   }
 594   return status;
 595 }
 596
 597 //______________________________________________________
 598 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
 599         //
 600         // Check ITS layer status
 601         //
 602         Int_t status = 0;
 603         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 604                 Int_t det;
 605                 Float_t xloc, zloc;
 606                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 607                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
 608         }
 609         return status;
 610 }
 611
 612
 613 //______________________________________________________
 614 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
 615   //
 616   // Checks if number of shared bits is above 1
 617   //
 618   Int_t nsharebit = 1;
 619   const TBits *shared;
 620   if(track->IsA() == AliESDtrack::Class())
 621     shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
 622   else if(track->IsA() == AliAODTrack::Class())
 623     shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
 624    else
 625     return 0;
 626
 627         if(shared->CountBits() >= 2) nsharebit=1;
 628   else nsharebit=0;
 629
 630   return nsharebit;
 631 }
 632
 633 //______________________________________________________
 634 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
 635         //
 636         // Get Number of findable clusters in the TPC
 637         //
 638         Int_t nClusters = 0; // in case no Information available consider all clusters findable
 639   TClass *type = track->IsA();
 640   if(TESTBIT(fTPCclusterDef, kFoundIter1)){
 641                 if(type == AliESDtrack::Class()){
 642       AliDebug(2, ("Using def kFoundIter1"));
 643       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 644       nClusters = esdtrack->GetTPCNclsIter1();
 645     } else {
 646       AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
 647     }
 648   } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
 649     AliDebug(2, ("Using def kCrossedRows"));
 650     nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
 651   } else{
 652     AliDebug(2, ("Using def kFound"));
 653           nClusters = track->GetTPCNcls();
 654   }
 655         return nClusters;
 656 }
 657
 658 //______________________________________________________
 659 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
 660   //
 661   // Get fraction of shared TPC clusters
 662   //
 663   Float_t fracClustersTPCShared = 0.0;
 664   Int_t nClustersTPC = track->GetTPCNcls();
 665   Int_t nClustersTPCShared = 0;
 666   TClass *type = track->IsA();
 667   if(type == AliESDtrack::Class()){
 668     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 669     nClustersTPCShared = esdtrack->GetTPCnclsS();
 670   } else if(type == AliESDtrack::Class()){
 671     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 672     const TBits &shared = aodtrack->GetTPCSharedMap();
 673     nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
 674   }
 675   if (nClustersTPC!=0) {
 676            fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
 677   }
 678   return fracClustersTPCShared;
 679 }
 680
 681 //______________________________________________________
 682 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
 683   //
 684   // Get Ratio of found / findable clusters for different definitions
 685   // Only implemented for ESD tracks
 686   //
 687   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
 688   TClass *type = track->IsA();
 689   if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
 690     AliDebug(2, "Using ratio def kCROverFindable");
 691     clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
 692   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
 693     AliDebug(2, "Using ratio def kFoundOverCR");
 694     clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
 695   } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
 696     if(type == AliESDtrack::Class()){
 697       AliDebug(2, "Using ratio def kFoundOverFindableIter1");
 698       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 699       clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
 700     } else {
 701       AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
 702       clusterRatio = 1.;
 703     }
 704   } else {
 705     AliDebug(2, "Using ratio def kFoundOverFindable");
 706     clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
 707   }
 708   return clusterRatio;
 709 }
 710
 711 //______________________________________________________
 712 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
 713   //
 714   // Get impact parameter
 715   //
 716   TClass *type = track->IsA();
 717   if(type == AliESDtrack::Class()){
 718     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 719     esdtrack->GetImpactParameters(radial, z);
 720   }
 721   else if(type == AliAODTrack::Class()){
 722     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 723     if(aodtrack){
 724       Double_t xyz[3];
 725       aodtrack->XYZAtDCA(xyz);
 726       //printf("xyz[0] %f, xyz[1] %f, xyz[2] %f\n",xyz[0], xyz[1],xyz[2]);
 727       if(xyz[0]==-999. || xyz[1]==-999. ||  xyz[2]==-999.){
 728         if(!fEvent) {
 729           z = xyz[2];
 730           radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
 731           //printf("No event\n");
 732         }
 733         else {
 734           // Propagate the global track to the DCA.
 735           const AliVVertex *vAOD = fEvent->GetPrimaryVertex();
 736           Double_t PosAtDCA[2] = {-999,-999};
 737           Double_t covar[3] = {-999,-999,-999};
 738           AliAODTrack *copyaodtrack = new AliAODTrack(*aodtrack);
 739           if(copyaodtrack->PropagateToDCA(vAOD,fEvent->GetMagneticField(),100.,PosAtDCA,covar)) {
 740             z = PosAtDCA[1];
 741             radial = PosAtDCA[0];
 742             //printf("Propagate z %f and radial %f\n",z,radial);
 743           }
 744           if(copyaodtrack) delete copyaodtrack;
 745         }
 746       }
 747       else {
 748         z = xyz[2];
 749         radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
 750       }
 751     }
 752   }
 753 }
 754 //______________________________________________________
 755 Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
 756   //
 757   // Is kink Daughter
 758   //
 759   TClass *type = track->IsA();
 760   if(type == AliESDtrack::Class()){
 761     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 762     if(!esdtrack) return kFALSE;
 763     if(esdtrack->GetKinkIndex(0)<=0) return kFALSE;
 764     else return kTRUE;
 765
 766   }
 767   else if(type == AliAODTrack::Class()){
 768     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 769     if(aodtrack){
 770       AliAODVertex *aodvertex = aodtrack->GetProdVertex();
 771       if(!aodvertex) return kFALSE;
 772       if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
 773       else return kFALSE;
 774     }
 775     else return kFALSE;
 776   }
 777
 778   return kFALSE;
 779 }
 780 //______________________________________________________
 781 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
 782   //
 783   // Is kink Mother: only for ESD since need to loop over vertices for AOD
 784   //
 785   //
 786
 787   TClass *type = track->IsA();
 788   if(type == AliESDtrack::Class()){
 789     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 790     if(!esdtrack) return kFALSE;
 791     if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
 792     else return kFALSE;
 793   }
 794
 795   return kFALSE;
 796
 797 }
 798 //______________________________________________________
 799 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
 800   //
 801   // Get ITS nb of clusters
 802   //
 803   TClass *type = track->IsA();
 804   if(type == AliESDtrack::Class()){
 805     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 806     return esdtrack->GetITSclusters(0);
 807
 808   }
 809   else if(type == AliAODTrack::Class()){
 810     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 811     if(aodtrack){
 812       return  aodtrack->GetITSNcls();
 813     }
 814   }
 815
 816   return 0;
 817 }
 818 //______________________________________________________
 819 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
 820         //
 821         // Get HFE impact parameter (with recalculated primary vertex)
 822         //
 823         dcaxy=0;
 824         dcansigmaxy=0;
 825   if(!fEvent){
 826     AliDebug(1, "No Input event available\n");
 827     return;
 828   }
 829   const Double_t kBeampiperadius=3.;
 830   TString type = track->IsA()->GetName();
 831   Double_t dcaD[2]={-999.,-999.};
 832   Double_t covD[3]={-999.,-999.,-999.};
 833   Float_t dcaF[2]={-999.,-999.};
 834   Float_t covF[3]={-999.,-999.,-999.};
 835
 836   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 837   if(!esdevent) {
 838     AliDebug(1, "No esd event available\n");
 839     return;
 840   }
 841   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
 842   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
 843    // recalculate primary vertex for peri. and pp
 844    AliVertexerTracks vertexer(fEvent->GetMagneticField());
 845    vertexer.SetITSMode();
 846    vertexer.SetMinClusters(4);
 847          Int_t skipped[2];
 848    skipped[0] = track->GetID();
 849    vertexer.SetSkipTracks(1,skipped);
 850
 851  //----diamond constraint-----------------------------
 852    vertexer.SetConstraintOn();
 853    Float_t diamondcovxy[3];
 854    esdevent->GetDiamondCovXY(diamondcovxy);
 855    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
 856    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 857    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 858    vertexer.SetVtxStart(diamond);
 859    delete diamond; diamond=NULL;
 860  //----------------------------------------------------
 861
 862    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 863
 864    // Getting the DCA
 865    // Propagation always done on a working copy to not disturb the track params of the original track
 866    AliESDtrack *esdtrack = NULL;
 867    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 868      // Case ESD track: take copy constructor
 869      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 870      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 871    } else {
 872     // Case AOD track: take different constructor
 873     esdtrack = new AliESDtrack(track);
 874    }
 875    if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
 876     // protection
 877     dcaxy = dcaD[0];
 878     if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
 879     if(!covD[0]) dcansigmaxy = -49.;
 880    }
 881    delete esdtrack;
 882    delete vtxESDSkip;
 883   }
 884   else{
 885    AliESDtrack *esdtrack = NULL;
 886    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 887     // Case ESD track: take copy constructor
 888     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 889     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 890    } else {
 891     // Case AOD track: take different constructor
 892     esdtrack = new AliESDtrack(track);
 893    }
 894    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
 895    dcaxy = dcaF[0];
 896    if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
 897    if(!covF[0]) dcansigmaxy = -49.;
 898    delete esdtrack;
 899   }
 900 }
 901
 902 //______________________________________________________
 903 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
 904         //
 905         // Get HFE impact parameter (with recalculated primary vertex)
 906         //
 907   if(!fEvent){
 908     AliDebug(1, "No Input event available\n");
 909     return;
 910   }
 911   const Double_t kBeampiperadius=3.;
 912   TString type = track->IsA()->GetName();
 913   Float_t dcaF[2]={-999.,-999.};
 914   Float_t covF[3]={-999.,-999.,-999.};
 915
 916   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 917   if(!esdevent) {
 918     AliDebug(1, "No esd event available\n");
 919     return;
 920   }
 921   const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
 922   if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
 923    // recalculate primary vertex for peri. and pp
 924    AliVertexerTracks vertexer(fEvent->GetMagneticField());
 925    vertexer.SetITSMode();
 926    vertexer.SetMinClusters(4);
 927          Int_t skipped[2];
 928    skipped[0] = track->GetID();
 929    vertexer.SetSkipTracks(1,skipped);
 930
 931  //----diamond constraint-----------------------------
 932    vertexer.SetConstraintOn();
 933    Float_t diamondcovxy[3];
 934    esdevent->GetDiamondCovXY(diamondcovxy);
 935    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
 936    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 937    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 938    vertexer.SetVtxStart(diamond);
 939    delete diamond; diamond=NULL;
 940  //----------------------------------------------------
 941
 942    vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 943
 944    // Getting the DCA
 945    // Propagation always done on a working copy to not disturb the track params of the original track
 946    AliESDtrack *esdtrack = NULL;
 947    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 948      // Case ESD track: take copy constructor
 949      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 950      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 951    } else {
 952     // Case AOD track: take different constructor
 953     esdtrack = new AliESDtrack(track);
 954    }
 955    if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
 956    delete esdtrack;
 957    delete vtxESDSkip;
 958   }
 959   else{
 960    AliESDtrack *esdtrack = NULL;
 961    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 962     // Case ESD track: take copy constructor
 963     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 964     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 965    } else {
 966     // Case AOD track: take different constructor
 967     esdtrack = new AliESDtrack(track);
 968    }
 969    if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
 970    dcaD[0] = dcaF[0];
 971    dcaD[1] = dcaF[1];
 972    covD[0] = covF[0];
 973    covD[1] = covF[1];
 974    covD[2] = covF[2];
 975    delete esdtrack;
 976   }
 977 }
 978
 979 //______________________________________________________
 980 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
 981         //
 982         // Get HFE impact parameter cut(pt dependent)
 983         //
 984
 985   TString type = track->IsA()->GetName();
 986   if(!type.CompareTo("AliESDtrack")){
 987     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 988     if(!esdtrack) return;
 989     Double_t pt = esdtrack->Pt();
 990     hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
 991     hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
 992   }
 993 }
 994 //______________________________________________________
 995 void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
 996         //
 997         // Get max impact parameter cut r (pt dependent)
 998         //
 999
1000   TString type = track->IsA()->GetName();
1001   if(!type.CompareTo("AliESDtrack")){
1002     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
1003     if(!esdtrack) return;
1004     Double_t pt = esdtrack->Pt();
1005     if(pt > 0.15) {
1006       maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
1007     }
1008     else maximpactRcut = 9999999999.0;
1009   }
1010 }