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     if(fQAlist) 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       if(fQAlist) 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   TH1 *htmp = NULL;
 311   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * kNhistos)))) htmp->Fill(impactR);
 312   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * kNhistos)))) htmp->Fill(impactZ);
 313   // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
 314   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * kNhistos)))) htmp->Fill(nTPCf > 0. ? static_cast<Float_t>(nclsTPC)/static_cast<Float_t>(nTPCf) : 1.);
 315   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * kNhistos)))) htmp->Fill(GetTRDnTrackletsPID(track));
 316   if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * kNhistos)))) htmp->Fill(nclsTPC);
 317   UChar_t itsPixel = track->GetITSClusterMap();
 318   TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * kNhistos));
 319   //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
 320   if(pixelHist){
 321     Double_t firstEntry = 0.5;
 322     if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
 323       pixelHist->Fill(firstEntry + 3);
 324     else{
 325       if(itsPixel & BIT(0)){
 326         pixelHist->Fill(firstEntry);
 327         if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
 328         else pixelHist->Fill(firstEntry + 4);
 329       }
 330       if(itsPixel & BIT(1)){
 331         pixelHist->Fill(firstEntry + 1);
 332         if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
 333       }
 334     }
 335   }
 336 }
 337
 338 // //______________________________________________________
 339 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
 340 //   //
 341 //   // Fill the QA histograms for MC tracks
 342 //   // Function can be called before cuts or after cut application (second argument)
 343 //   // Not yet implemented
 344 //   //
 345 // }
 346
 347 //______________________________________________________
 348 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
 349   //
 350   // Fill cut correlation histograms for tracks that didn't pass cuts
 351   //
 352   const Int_t kNhistos = 6;
 353   TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * kNhistos));
 354   if(correlation){
 355     for(Int_t icut = 0; icut < kNcuts; icut++){
 356       if(!TESTBIT(fRequirements, icut)) continue;
 357       for(Int_t jcut = icut; jcut < kNcuts; jcut++){
 358         if(!TESTBIT(fRequirements, jcut)) continue;
 359         if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
 360                 correlation->Fill(icut, jcut);
 361       }
 362     }
 363   }
 364 }
 365
 366 //______________________________________________________
 367 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
 368   //
 369   // Add QA histograms
 370   // For each cut a histogram before and after track cut is created
 371   // Histos before respectively after cut are stored in different lists
 372   // Additionally a histogram with the cut correlation is created and stored
 373   // in the top directory
 374   //
 375
 376   const Int_t kNhistos = 6;
 377   TH1 *histo1D = 0x0;
 378   TH2 *histo2D = 0x0;
 379   TString cutstr[2] = {"before", "after"};
 380
 381   if(!fQAlist) fQAlist = new TList;  // for internal representation, not owner
 382   for(Int_t icond = 0; icond < 2; icond++){
 383     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * kNhistos);
 384     fQAlist->AddAt(histo1D, 0 + icond * kNhistos);
 385     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 386     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 387     qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * kNhistos);
 388     fQAlist->AddAt(histo1D, 1 + icond * kNhistos);
 389     histo1D->GetXaxis()->SetTitle("Impact Parameter");
 390     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 391     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2 + icond * kNhistos);
 392     fQAlist->AddAt(histo1D, 2 + icond * kNhistos);
 393     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
 394     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 395     qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * kNhistos);
 396     fQAlist->AddAt(histo1D, 3 + icond * kNhistos);
 397     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
 398     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 399     qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * kNhistos);
 400     fQAlist->AddAt(histo1D, 4 + icond * kNhistos);
 401     histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
 402     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 403     qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * kNhistos);
 404     fQAlist->AddAt(histo1D, 5 + icond * kNhistos);
 405     histo1D->GetXaxis()->SetTitle("ITS Pixel");
 406     histo1D->GetYaxis()->SetTitle("Number of Tracks");
 407     Int_t first = histo1D->GetXaxis()->GetFirst();
 408     TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
 409     for(Int_t ilabel = 0; ilabel < 6; ilabel++)
 410       histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
 411   }
 412   // Add cut correlation
 413   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * kNhistos);
 414   fQAlist->AddAt(histo2D, 2 * kNhistos);
 415   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters"};
 416   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
 417   for(Int_t icut = 0; icut < kNcuts; icut++){
 418     histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
 419     histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
 420   }
 421 }
 422
 423 //______________________________________________________
 424 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
 425   for(Int_t ibit = 32; ibit--; )
 426     printf("%d", bitmap & BIT(ibit) ? 1 : 0);
 427   printf("\n");
 428 }
 429
 430 //______________________________________________________
 431 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 432   //
 433   // Check whether ITS area is dead
 434   //
 435   Bool_t status;
 436   switch(itsStatus){
 437     case 2: status = kFALSE; break;
 438     case 3: status = kFALSE; break;
 439     case 7: status = kFALSE; break;
 440     default: status = kTRUE;
 441   }
 442   return status;
 443 }
 444
 445 //______________________________________________________
 446 Int_t AliHFEextraCuts::GetTRDnTrackletsPID(AliVTrack *track){
 447         //
 448         // Get Number of TRD tracklets
 449         //
 450         Int_t nTracklets = 0;
 451         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 452                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 453                 if(esdtrack) nTracklets = esdtrack->GetTRDntrackletsPID();
 454         } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODTrack")){
 455                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 456                 AliAODPid *pidobject = NULL;
 457     if(aodtrack) pidobject = aodtrack->GetDetPid();
 458                 // this is normally NOT the way to do this, but due to limitation in the
 459                 // AOD track it is not possible in a different way
 460                 if(pidobject){
 461                         Float_t *trdmom = pidobject->GetTRDmomentum();
 462                         for(Int_t ily = 0; ily < 6; ily++){
 463                                 if(trdmom[ily] > -1) nTracklets++;
 464                         }
 465                 } else nTracklets = 6;  // No Cut possible
 466         }
 467         return nTracklets;
 468 }
 469
 470 //______________________________________________________
 471 Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
 472         //
 473         // Check ITS layer status
 474         //
 475         Int_t status = 0;
 476         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 477                 Int_t det;
 478                 Float_t xloc, zloc;
 479                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 480                 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
 481         }
 482         return status;
 483 }
 484
 485 //______________________________________________________
 486 Int_t AliHFEextraCuts::GetTPCfindableClusters(AliVTrack *track, Bool_t iter1){
 487         //
 488         // Get Number of findable clusters in the TPC
 489         //
 490   AliDebug(1, Form("Using TPC clusters from iteration 1: %s", iter1 ? "Yes" : "No"));
 491         Int_t nClusters = 159; // in case no Information available consider all clusters findable
 492         if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 493     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 494                 if(esdtrack) nClusters = esdtrack->GetTPCNclsF();
 495   }
 496         return nClusters;
 497 }
 498
 499 //______________________________________________________
 500 Int_t AliHFEextraCuts::GetTPCncls(AliVTrack *track, Bool_t iter1){
 501         //
 502         // Get Number of findable clusters in the TPC
 503         //
 504   AliDebug(1, Form("Using TPC clusters from iteration 1: %s", iter1 ? "Yes" : "No"));
 505         Int_t nClusters = 0; // in case no Information available consider all clusters findable
 506         TString type = track->IsA()->GetName();
 507         if(!type.CompareTo("AliESDtrack")){
 508     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 509     if(esdtrack){ // coverity
 510       if(iter1)
 511                     nClusters = esdtrack->GetTPCNclsIter1();
 512       else
 513                     nClusters = esdtrack->GetTPCNcls();
 514     }
 515   }
 516         else if(!type.CompareTo("AliAODTrack")){
 517                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 518           const TBits &tpcmap = aodtrack->GetTPCClusterMap();
 519           for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
 520                   if(tpcmap.TestBitNumber(ibit)) nClusters++;
 521         }
 522         return nClusters;
 523 }
 524
 525 //______________________________________________________
 526 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
 527         //
 528         // Get impact parameter
 529         //
 530         TString type = track->IsA()->GetName();
 531         if(!type.CompareTo("AliESDtrack")){
 532                 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 533     if(esdtrack) esdtrack->GetImpactParameters(radial, z);
 534         }
 535         else if(!type.CompareTo("AliAODTrack")){
 536                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 537     if(aodtrack){
 538                   Double_t xyz[3];
 539                   aodtrack->XYZAtDCA(xyz);
 540                   z = xyz[2];
 541                   radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
 542     }
 543         }
 544 }
 545
 546 //______________________________________________________
 547 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
 548         //
 549         // Get HFE impact parameter (with recalculated primary vertex)
 550         //
 551         dcaxy=0;
 552         dcansigmaxy=0;
 553   if(!fEvent){
 554     AliDebug(1, "No Input event available\n");
 555     return;
 556   }
 557   const Double_t kBeampiperadius=3.;
 558   TString type = track->IsA()->GetName();
 559   Double_t dca[2]={-999.,-999.};
 560   Double_t cov[3]={-999.,-999.,-999.};
 561
 562   // recalculate primary vertex
 563   AliVertexerTracks vertexer(fEvent->GetMagneticField());
 564   vertexer.SetITSMode();
 565   vertexer.SetMinClusters(4);
 566         Int_t skipped[2];
 567   skipped[0] = track->GetID();
 568   vertexer.SetSkipTracks(1,skipped);
 569   AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 570   vertexer.SetSkipTracks(1,skipped);
 571   if(vtxESDSkip->GetNContributors()<2) return;
 572
 573   // Getting the DCA
 574   // Propagation always done on a working copy to not disturb the track params of the original track
 575   AliESDtrack *esdtrack = NULL;
 576   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 577     // Case ESD track: take copy constructor
 578     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 579     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 580   } else {
 581     // Case AOD track: take different constructor
 582     esdtrack = new AliESDtrack(track);
 583   }
 584   if(esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dca, cov)){
 585     // protection
 586     dcaxy = dca[0];
 587     if(cov[0]) dcansigmaxy = dcaxy/TMath::Sqrt(cov[0]);
 588     if(!cov[0]) dcansigmaxy = -99.;
 589   }
 590   delete esdtrack;
 591   delete vtxESDSkip;
 592 }
 593
 594
 595 //______________________________________________________
 596 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
 597         //
 598         // Get HFE impact parameter cut(pt dependent)
 599         //
 600
 601   TString type = track->IsA()->GetName();
 602   if(!type.CompareTo("AliESDtrack")){
 603     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 604     if(!esdtrack) return;
 605     Double_t pt = esdtrack->Pt();
 606     //hfeimpactRcut=0.0064+0.078*exp(-0.56*pt);  // used Carlo's old parameter
 607     hfeimpactRcut=0.011+0.077*exp(-0.65*pt); // used Carlo's new parameter
 608     hfeimpactnsigmaRcut=3; // 3 sigma trail cut
 609   }
 610 }
 611