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