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