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