PWG3/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 "AliESDEvent.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDtrack.h"
  38 #include "AliLog.h"
  39 #include "AliMCParticle.h"
  40 #include "AliVEvent.h"
  41 #include "AliVTrack.h"
  42 #include "AliVParticle.h"
  43 #include "AliVertexerTracks.h"
  44 #include "AliVVertex.h"
  45
  46 #include "AliHFEextraCuts.h"
  47
  48 ClassImp(AliHFEextraCuts)
  49
  50 //______________________________________________________
  51 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
  52   AliCFCutBase(name, title),
  53   fEvent(NULL),
  54   fCutCorrelation(0),
  55   fRequirements(0),
  56   fMinNClustersTPC(0),
  57   fClusterRatioTPC(0.),
  58   fMinTrackletsTRD(0),
  59   fPixelITS(0),
  60   fTOFpid(kFALSE),
  61   fTOFmismatch(kFALSE),
  62   fTPCclusterDef(0),
  63   fTPCclusterRatioDef(0),
  64   fCheck(kFALSE),
  65   fQAlist(0x0) ,
  66   fDebugLevel(0)
  67 {
  68   //
  69   // Default Constructor
  70   //
  71   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
  72   memset(fIPcutParam, 0, sizeof(Float_t) * 4);
  73 }
  74
  75 //______________________________________________________
  76 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
  77   AliCFCutBase(c),
  78   fEvent(c.fEvent),
  79   fCutCorrelation(c.fCutCorrelation),
  80   fRequirements(c.fRequirements),
  81   fMinNClustersTPC(c.fMinNClustersTPC),
  82   fClusterRatioTPC(c.fClusterRatioTPC),
  83   fMinTrackletsTRD(c.fMinTrackletsTRD),
  84   fPixelITS(c.fPixelITS),
  85   fTOFpid(c.fTOFpid),
  86   fTOFmismatch(c.fTOFmismatch),
  87   fTPCclusterDef(c.fTPCclusterDef),
  88   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
  89   fCheck(c.fCheck),
  90   fQAlist(0x0),
  91   fDebugLevel(0)
  92 {
  93   //
  94   // Copy constructor
  95   // Performs a deep copy
  96   //
  97   memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
  98   memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
  99   if(IsQAOn()){
 100     fIsQAOn = kTRUE;
 101     fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 102     if(fQAlist) fQAlist->SetOwner();
 103   }
 104 }
 105
 106 //______________________________________________________
 107 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
 108   //
 109   // Assignment operator
 110   //
 111   if(this != &c){
 112     AliCFCutBase::operator=(c);
 113     fEvent = c.fEvent;
 114     fCutCorrelation = c.fCutCorrelation;
 115     fRequirements = c.fRequirements;
 116     fClusterRatioTPC = c.fClusterRatioTPC;
 117     fMinNClustersTPC = c.fMinNClustersTPC;
 118     fMinTrackletsTRD = c.fMinTrackletsTRD;
 119     fPixelITS = c.fPixelITS;
 120     fTPCclusterDef = c.fTPCclusterDef;
 121     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
 122     fTOFpid = c.fTOFpid;
 123     fTOFmismatch = c.fTOFmismatch;
 124     fCheck = c.fCheck;
 125     fDebugLevel = c.fDebugLevel;
 126     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
 127     memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
 128     if(IsQAOn()){
 129       fIsQAOn = kTRUE;
 130       fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
 131       if(fQAlist) fQAlist->SetOwner();
 132     }else fQAlist = 0x0;
 133   }
 134   return *this;
 135 }
 136
 137 //______________________________________________________
 138 AliHFEextraCuts::~AliHFEextraCuts(){
 139   //
 140   // Destructor
 141   //
 142   if(fQAlist) delete fQAlist;
 143 }
 144
 145 //______________________________________________________
 146 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
 147   //
 148   // Set Virtual event an make a copy
 149   //
 150   if (!event) {
 151     AliError("Pointer to AliVEvent !");
 152     return;
 153   }
 154   TString className(event->ClassName());
 155   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
 156     AliError("argument must point to an AliESDEvent or AliAODEvent !");
 157     return ;
 158   }
 159   fEvent = (AliVEvent*) event;
 160
 161 }
 162
 163 //______________________________________________________
 164 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
 165   //
 166   // Steering function for the track selection
 167   //
 168   TString type = o->IsA()->GetName();
 169   AliDebug(2, Form("Object type %s", type.Data()));
 170   if(!type.CompareTo("AliESDtrack") || !type.CompareTo("AliAODTrack")){
 171     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
 172   }
 173   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
 174 }
 175
 176 //______________________________________________________
 177 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
 178   //
 179   // Checks cuts on reconstructed tracks
 180   // returns true if track is selected
 181   // QA histograms are filled before track selection and for
 182   // selected tracks after track selection
 183   //
 184   AliDebug(1, "Called");
 185   ULong64_t survivedCut = 0;    // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
 186   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
 187   // Apply cuts
 188   Float_t impactR, impactZ;
 189   Double_t hfeimpactR, hfeimpactnsigmaR;
 190   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
 191   Bool_t tofstep = kTRUE;
 192   Bool_t tofmismatchstep = kTRUE;
 193   GetImpactParameters(track, impactR, impactZ);
 194   if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 195     // Protection for PbPb
 196     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
 197     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
 198   }
 199   UInt_t nclsTPC = GetTPCncls(track);
 200   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 201   Double_t ratioTPC = GetTPCclusterRatio(track);
 202   UChar_t trdTracklets;
 203   trdTracklets = GetTRDnTrackletsPID(track);
 204   UChar_t itsPixel = track->GetITSClusterMap();
 205   Int_t status1 = GetITSstatus(track, 0);
 206   Int_t status2 = GetITSstatus(track, 1);
 207   Bool_t statusL0 = CheckITSstatus(status1);
 208   Bool_t statusL1 = CheckITSstatus(status2);
 209   if(TESTBIT(fRequirements, kMinImpactParamR)){
 210     // cut on min. Impact Parameter in Radial direction
 211     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
 212   }
 213   if(TESTBIT(fRequirements, kMinImpactParamZ)){
 214     // cut on min. Impact Parameter in Z direction
 215     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
 216   }
 217   if(TESTBIT(fRequirements, kMaxImpactParamR)){
 218     // cut on max. Impact Parameter in Radial direction
 219     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
 220   }
 221   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
 222     // cut on max. Impact Parameter in Z direction
 223     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
 224   }
 225   if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
 226     // cut on min. HFE Impact Parameter in Radial direction
 227     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
 228   }
 229   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
 230     // cut on max. HFE Impact Parameter n sigma in Radial direction
 231     if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
 232   }
 233   if(TESTBIT(fRequirements, kClusterRatioTPC)){
 234     // cut on min ratio of found TPC clusters vs findable TPC clusters
 235     if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
 236   }
 237   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
 238     // cut on minimum number of TRD tracklets
 239     AliDebug(1, Form("Min TRD cut: [%d|%d]\n", fMinTrackletsTRD, trdTracklets));
 240     if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
 241   }
 242   if(TESTBIT(fRequirements, kMinNClustersTPC)){
 243     // cut on minimum number of TRD tracklets
 244     AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
 245     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
 246   }
 247   if(TESTBIT(fRequirements, kPixelITS)){
 248     // cut on ITS pixel layers
 249     AliDebug(1, "ITS cluster Map: ");
 250     //PrintBitMap(itsPixel);
 251     switch(fPixelITS){
 252       case kFirst:
 253         AliDebug(2, "First");
 254               if(TESTBIT(itsPixel, 0))
 255                 SETBIT(survivedCut, kPixelITS);
 256         else
 257                 if(fCheck && !statusL0)
 258                   SETBIT(survivedCut, kPixelITS);
 259                     break;
 260       case kSecond:
 261         AliDebug(2, "Second");
 262               if(TESTBIT(itsPixel, 1))
 263                 SETBIT(survivedCut, kPixelITS);
 264         else
 265                 if(fCheck && !statusL1)
 266                   SETBIT(survivedCut, kPixelITS);
 267                     break;
 268       case kBoth:
 269         AliDebug(2, "Both");
 270               if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
 271                       SETBIT(survivedCut, kPixelITS);
 272               else
 273           if(fCheck && !(statusL0 && statusL1))
 274                         SETBIT(survivedCut, kPixelITS);
 275               break;
 276       case kAny:
 277         AliDebug(2, "Any");
 278               if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
 279                 SETBIT(survivedCut, kPixelITS);
 280         else
 281                 if(fCheck && !(statusL0 || statusL1))
 282                     SETBIT(survivedCut, kPixelITS);
 283                     break;
 284       default:
 285         AliDebug(2, "None");
 286         break;
 287     }
 288     AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
 289   }
 290   if(fTOFpid){
 291     // cut on TOF pid
 292     tofstep = kFALSE;
 293     if(track->GetStatus() & AliESDtrack::kTOFpid) tofstep = kTRUE;
 294
 295   }
 296
 297   if(fTOFmismatch){
 298     // cut on TOF mismatch
 299     tofmismatchstep = kFALSE;
 300     if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) tofmismatchstep = kTRUE;
 301
 302   }
 303
 304
 305   if(fRequirements == survivedCut && tofstep && tofmismatchstep){
 306     //
 307     // Track selected
 308     //
 309     AliDebug(2, "Track Survived cuts\n");
 310     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
 311     return kTRUE;
 312   }
 313   AliDebug(2, "Track cut");
 314   if(IsQAOn()) FillCutCorrelation(survivedCut);
 315   return kFALSE;
 316 }
 317
 318 //______________________________________________________
 319 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
 320   //
 321   // Checks cuts on Monte Carlo tracks
 322   // returns true if track is selected
 323   // QA histograms are filled before track selection and for
 324   // selected tracks after track selection
 325   //
 326   return kTRUE; // not yet implemented
 327 }
 328
 329 //______________________________________________________
 330 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
 331   //
 332   // Fill the QA histograms for ESD tracks
 333   // Function can be called before cuts or after cut application (second argument)
 334   //
 335   const Int_t kNhistos = 6;
 336   Float_t impactR, impactZ;
 337   GetImpactParameters(track, impactR, impactZ);
 338   TH1 *htmp = NULL;
 339   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * kNhistos)))) htmp->Fill(impactR);
 340   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * kNhistos)))) htmp->Fill(impactZ);
 341   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 342   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * kNhistos)))) htmp->Fill(GetTPCclusterRatio(track));
 343   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * kNhistos)))) htmp->Fill(GetTRDnTrackletsPID(track));
 344   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * kNhistos)))) htmp->Fill(GetTPCncls(track));
 345   UChar_t itsPixel = track->GetITSClusterMap();
 346   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * kNhistos));
 347   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 348   if(pixelHist){
 349     Double_t firstEntry = 0.5;
 350     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 351       pixelHist->Fill(firstEntry + 3);
 352     else{
 353       if(itsPixel & BIT(0)){
 354         pixelHist->Fill(firstEntry);
 355         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 356         else pixelHist->Fill(firstEntry + 4);
 357       }
 358       if(itsPixel & BIT(1)){
 359         pixelHist->Fill(firstEntry + 1);
 360         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 361       }
 362     }
 363   }
 364 }
 365
 366 // //______________________________________________________
 367 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 368 //   //
 369 //   // Fill the QA histograms for MC tracks
 370 //   // Function can be called before cuts or after cut application (second argument)
 371 //   // Not yet implemented
 372 //   //
 373 // }
 374
 375 //______________________________________________________
 376 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 377   //
 378   // Fill cut correlation histograms for tracks that didn't pass cuts
 379   //
 380   const Int_t kNhistos = 6;
 381   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * kNhistos));
 382   if(correlation){
 383     for(Int_t icut = 0; icut < kNcuts; icut++){
 384       if(!TESTBIT(fRequirements, icut)) continue;
 385       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 386         if(!TESTBIT(fRequirements, jcut)) continue;
 387         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 388                 correlation->Fill(icut, jcut);
 389       }
 390     }
 391   }
 392 }
 393
 394 //______________________________________________________
 395 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 396   //
 397   // Add QA histograms
 398   // For each cut a histogram before and after track cut is created
 399   // Histos before respectively after cut are stored in different lists
 400   // Additionally a histogram with the cut correlation is created and stored
 401   // in the top directory
 402   //
 403
 404   const Int_t kNhistos = 6;
 405   TH1 *histo1D = 0x0;
 406   TH2 *histo2D = 0x0;
 407   TString cutstr[2] = {"before", "after"};
 408
 409   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
 410   for(Int_t icond = 0; icond < 2; icond++){
 411     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * kNhistos);
 412     fQAlist->AddAt(histo1D, 0 + icond * kNhistos);
 413     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 414     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 415     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * kNhistos);
 416     fQAlist->AddAt(histo1D, 1 + icond * kNhistos);
 417     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 418     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 419     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * kNhistos);
 420     fQAlist->AddAt(histo1D, 2 + icond * kNhistos);
 421     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 422     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 423     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * kNhistos);
 424     fQAlist->AddAt(histo1D, 3 + icond * kNhistos);
 425     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 426     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 427     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * kNhistos);
 428     fQAlist->AddAt(histo1D, 4 + icond * kNhistos);
 429     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
 430     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 431     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * kNhistos);
 432     fQAlist->AddAt(histo1D, 5 + icond * kNhistos);
 433     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 434     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 435     Int_t first = histo1D->GetXaxis()->GetFirst();
 436     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 437     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 438       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 439   }
 440   // Add cut correlation
 441   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * kNhistos);
 442   fQAlist->AddAt(histo2D, 2 * kNhistos);
 443   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters"};
 444   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 445   for(Int_t icut = 0; icut < kNcuts; icut++){
 446     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 447     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 448   }
 449 }
 450
 451 //______________________________________________________
 452 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 453   for(Int_t ibit = 32; ibit--; )
 454     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 455   printf("\n");
 456 }
 457
 458 //______________________________________________________
 459 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 460   //
 461   // Check whether ITS area is dead
 462   //
 463   Bool_t status;
 464   switch(itsStatus){
 465     case 2: status = kFALSE; break;
 466     case 3: status = kFALSE; break;
 467     case 7: status = kFALSE; break;
 468     default: status = kTRUE;
 469   }
 470   return status;
 471 }
 472
 473 //______________________________________________________
 474 Int_t AliHFEextraCuts::GetTRDnTrackletsPID(AliVTrack *track){
 475         //
 476         // Get Number of TRD tracklets
 477         //
 478         Int_t nTracklets = 0;
 479         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 480                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 481                 if(esdtrack) nTracklets = esdtrack->GetTRDntrackletsPID();
 482         } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODTrack")){
 483                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 484                 AliAODPid *pidobject = NULL;
 485     if(aodtrack) pidobject = aodtrack->GetDetPid();
 486                 // this is normally NOT the way to do this, but due to limitation in the
 487                 // AOD track it is not possible in a different way
 488                 if(pidobject){
 489                         Float_t *trdmom = pidobject->GetTRDmomentum();
 490                         for(Int_t ily = 0; ily < 6; ily++){
 491                                 if(trdmom[ily] > -1) nTracklets++;
 492                         }
 493                 } else nTracklets = 6;  // No Cut possible
 494         }
 495         return nTracklets;
 496 }
 497
 498 //______________________________________________________
 499 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
 500         //
 501         // Check ITS layer status
 502         //
 503         Int_t status = 0;
 504         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 505                 Int_t det;
 506                 Float_t xloc, zloc;
 507                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 508                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
 509         }
 510         return status;
 511 }
 512
 513 //______________________________________________________
 514 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
 515         //
 516         // Get Number of findable clusters in the TPC
 517         //
 518         Int_t nClusters = 0; // in case no Information available consider all clusters findable
 519         TString type = track->IsA()->GetName();
 520         if(!type.CompareTo("AliESDtrack")){
 521     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 522     if(esdtrack){ // coverity
 523       if(TESTBIT(fTPCclusterDef, kFoundIter1)){
 524                     nClusters = esdtrack->GetTPCNclsIter1();
 525         AliDebug(2, ("Using def kFoundIter1"));
 526       } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
 527         AliDebug(2, ("Using def kCrossedRows"));
 528         nClusters = static_cast<UInt_t>(esdtrack->GetTPCClusterInfo(2,1));
 529       } else{
 530         AliDebug(2, ("Using def kFound"));
 531                     nClusters = esdtrack->GetTPCNcls();
 532       }
 533     }
 534   }
 535         else if(!type.CompareTo("AliAODTrack")){
 536                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 537     if(aodtrack){
 538             const TBits &tpcmap = aodtrack->GetTPCClusterMap();
 539             for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
 540                     if(tpcmap.TestBitNumber(ibit)) nClusters++;
 541     }
 542         }
 543         return nClusters;
 544 }
 545
 546 //______________________________________________________
 547 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
 548   //
 549   // Get Ratio of found / findable clusters for different definitions
 550   // Only implemented for ESD tracks
 551   //
 552   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
 553   TString type = track->IsA()->GetName();
 554   if(!type.CompareTo("AliESDtrack")){
 555     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 556     if(esdtrack){ // coverity
 557       if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
 558         AliDebug(2, "Using ratio def kCROverFindable");
 559         clusterRatio = esdtrack->GetTPCNclsF() ? esdtrack->GetTPCClusterInfo(2,1)/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // crossed rows/findable
 560       } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
 561         AliDebug(2, "Using ratio def kFoundOverCR");
 562         clusterRatio = esdtrack->GetTPCClusterInfo(2,0); // found/crossed rows
 563       } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
 564         AliDebug(2, "Using ratio def kFoundOverFindableIter1");
 565         clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
 566       } else {
 567         AliDebug(2, "Using ratio def kFoundOverFindable");
 568         clusterRatio = esdtrack->GetTPCNclsF() ? static_cast<Double_t>(esdtrack->GetTPCNcls())/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // found/findable
 569       }
 570     }
 571   }
 572   return clusterRatio;
 573 }
 574
 575 //______________________________________________________
 576 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
 577         //
 578         // Get impact parameter
 579         //
 580         TString type = track->IsA()->GetName();
 581         if(!type.CompareTo("AliESDtrack")){
 582                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 583     if(esdtrack) esdtrack->GetImpactParameters(radial, z);
 584         }
 585         else if(!type.CompareTo("AliAODTrack")){
 586                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 587     if(aodtrack){
 588                   Double_t xyz[3];
 589                   aodtrack->XYZAtDCA(xyz);
 590                   z = xyz[2];
 591                   radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
 592     }
 593         }
 594 }
 595
 596 //______________________________________________________
 597 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
 598         //
 599         // Get HFE impact parameter (with recalculated primary vertex)
 600         //
 601         dcaxy=0;
 602         dcansigmaxy=0;
 603   if(!fEvent){
 604     AliDebug(1, "No Input event available\n");
 605     return;
 606   }
 607   const Double_t kBeampiperadius=3.;
 608   TString type = track->IsA()->GetName();
 609   Double_t dca[2]={-999.,-999.};
 610   Double_t cov[3]={-999.,-999.,-999.};
 611
 612   // recalculate primary vertex
 613   AliVertexerTracks vertexer(fEvent->GetMagneticField());
 614   vertexer.SetITSMode();
 615   vertexer.SetMinClusters(4);
 616         Int_t skipped[2];
 617   skipped[0] = track->GetID();
 618   vertexer.SetSkipTracks(1,skipped);
 619   AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 620   vertexer.SetSkipTracks(1,skipped);
 621   if(vtxESDSkip->GetNContributors()<2) return;
 622
 623   // Getting the DCA
 624   // Propagation always done on a working copy to not disturb the track params of the original track
 625   AliESDtrack *esdtrack = NULL;
 626   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 627     // Case ESD track: take copy constructor
 628     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 629     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 630   } else {
 631     // Case AOD track: take different constructor
 632     esdtrack = new AliESDtrack(track);
 633   }
 634   if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dca, cov)){
 635     // protection
 636     dcaxy = dca[0];
 637     if(cov[0]) dcansigmaxy = dcaxy/TMath::Sqrt(cov[0]);
 638     if(!cov[0]) dcansigmaxy = -99.;
 639   }
 640   delete esdtrack;
 641   delete vtxESDSkip;
 642 }
 643
 644
 645 //______________________________________________________
 646 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
 647         //
 648         // Get HFE impact parameter cut(pt dependent)
 649         //
 650
 651   TString type = track->IsA()->GetName();
 652   if(!type.CompareTo("AliESDtrack")){
 653     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 654     if(!esdtrack) return;
 655     Double_t pt = esdtrack->Pt();
 656     hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
 657     hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
 658   }
 659 }