PWGHF/hfe/AliSelectNonHFE.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
  17 ////////////////////////////////////////////////////////////////////////
  18 //                                                                    //
  19 //  Class for the Selection of Non-Heavy-Flavour-Electrons trought        //
  20 //      the invariant mass method. The selection can be done from two     //
  21 //  diferent algorithms, which can be choosed calling the function    //
  22 //  "SetAlgorithm(TString Algorithm)".                                                            //
  23 //                                                                    //
  24 //              Author: Elienos Pereira de Oliveira Filho                                 //
  25 //                                      (University of São Paulo)                            //
  26 //                                                                    //
  27 ////////////////////////////////////////////////////////////////////////
  28
  29 #include "TH1F.h"
  30 #include "TMath.h"
  31 #include "TLorentzVector.h"
  32 #include "AliAODTrack.h"
  33 #include "AliESDtrack.h"
  34 #include "AliVEvent.h"
  35 #include "AliESDtrackCuts.h"
  36 #include "AliVTrack.h"
  37 #include "AliVParticle.h"
  38 #include "AliPIDResponse.h"
  39 #include "AliPID.h"
  40 #include "AliExternalTrackParam.h"
  41 #include "AliSelectNonHFE.h"
  42 #include "AliKFParticle.h"
  43 #include "AliLog.h"
  44 #include "stdio.h"
  45 #include "iostream"
  46 #include "fstream"
  47
  48 ClassImp(AliSelectNonHFE)
  49 //________________________________________________________________________
  50 AliSelectNonHFE::AliSelectNonHFE(const char *name, const Char_t *title)
  51 : TNamed(name, title)
  52   ,fTrackCuts(0)
  53   ,fAlgorithm("DCA")
  54   ,fAngleCut(999)
  55   ,fdcaCut(999)
  56   ,fTPCnSigmaMin(-3)
  57   ,fTPCnSigmaMax(3)
  58   ,fMassCut(0.5)
  59   ,fChi2OverNDFCut(999)
  60   ,fPtMin(0.3)
  61   ,fIsLS(kFALSE)
  62   ,fIsULS(kFALSE)
  63   ,fIsAOD(kFALSE)
  64   ,fHasPtCut(kFALSE)
  65   ,fNLS(0)
  66   ,fNULS(0)
  67   ,fTpcNcls(50)
  68   ,fLSPartner(0)
  69   ,fULSPartner(0)
  70   ,fHistMass(0)
  71   ,fHistMassBack(0)
  72   ,fHistDCA(0)
  73   ,fHistDCABack(0)
  74   ,fHistAngle(0)
  75   ,fHistAngleBack(0)
  76   ,fPIDResponse(0)
  77 {
  78   //
  79   // Constructor
  80   //
  81
  82   fTrackCuts = new AliESDtrackCuts();
  83   //Configure Default Track Cuts
  84   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  85   fTrackCuts->SetRequireTPCRefit(kTRUE);
  86   fTrackCuts->SetEtaRange(-0.9,0.9);
  87   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
  88   fTrackCuts->SetMaxChi2PerClusterTPC(4.0);
  89   fTrackCuts->SetMinNClustersTPC(50);
  90   fTrackCuts->SetPtRange(0.3,1e10);
  91
  92
  93 }
  94
  95 //________________________________________________________________________
  96 AliSelectNonHFE::AliSelectNonHFE()
  97   : TNamed()
  98   ,fTrackCuts(0)
  99   ,fAlgorithm("DCA")
 100   ,fAngleCut(999)
 101   ,fdcaCut(999)
 102   ,fTPCnSigmaMin(-3)
 103   ,fTPCnSigmaMax(3)
 104   ,fMassCut(0.5)
 105   ,fChi2OverNDFCut(999)
 106   ,fPtMin(0.3)
 107   ,fIsLS(kFALSE)
 108   ,fIsULS(kFALSE)
 109   ,fIsAOD(kFALSE)
 110   ,fHasPtCut(kFALSE)
 111   ,fNLS(0)
 112   ,fNULS(0)
 113   ,fTpcNcls(50)
 114   ,fLSPartner(0)
 115   ,fULSPartner(0)
 116   ,fHistMass(0)
 117   ,fHistMassBack(0)
 118   ,fHistDCA(0)
 119   ,fHistDCABack(0)
 120   ,fHistAngle(0)
 121   ,fHistAngleBack(0)
 122   ,fPIDResponse(0)
 123 {
 124   //
 125   // Constructor
 126   //
 127
 128   fTrackCuts = new AliESDtrackCuts();
 129   //Configure Default Track Cuts
 130   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 131   fTrackCuts->SetRequireTPCRefit(kTRUE);
 132   fTrackCuts->SetEtaRange(-0.9,0.9);
 133   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 134   fTrackCuts->SetMaxChi2PerClusterTPC(4.0);
 135   fTrackCuts->SetMinNClustersTPC(50);
 136   fTrackCuts->SetPtRange(0.3,1e10);
 137
 138
 139 }
 140
 141 //_________________________________________
 142 AliSelectNonHFE::~AliSelectNonHFE()
 143 {
 144   //
 145   // Destructor
 146   //
 147
 148   if(fTrackCuts) delete fTrackCuts;
 149   if(fLSPartner) delete [] fLSPartner;
 150   if(fULSPartner) delete [] fULSPartner;
 151
 152
 153 }
 154
 155 //__________________________________________
 156 void AliSelectNonHFE::FindNonHFE(Int_t iTrack1, AliVParticle *Vtrack1, AliVEvent *fVevent)
 157 {
 158   AliVTrack *track1 = dynamic_cast<AliVTrack*>(Vtrack1);
 159   AliESDtrack *etrack1 = dynamic_cast<AliESDtrack*>(Vtrack1);
 160   AliExternalTrackParam *extTrackParam1=NULL;
 161   AliExternalTrackParam *extTrackParam2=NULL;
 162   if(fAlgorithm=="DCA")
 163   {
 164         extTrackParam1 = new AliExternalTrackParam();
 165     extTrackParam1->CopyFromVTrack(track1);
 166         }
 167   //
 168   // Find non HFE electrons
 169   //
 170
 171   //Magnetic Field
 172   Double_t bfield = fVevent->GetMagneticField();
 173
 174   //Second Track loop
 175
 176   fIsULS = kFALSE;      //Non-HFE Unlike signal Flag
 177   fIsLS = kFALSE;       //Non-HFE like signal Flag
 178   fNULS = 0;            //Non-HFE Unlike signal Flag
 179   fNLS = 0;             //Non-HFE like signal Flag
 180
 181   if(fLSPartner) delete [] fLSPartner;
 182   if(fULSPartner) delete [] fULSPartner;
 183   fLSPartner = new int [100];   //store the partners index
 184   fULSPartner = new int [100];  //store the partners index
 185
 186   for(Int_t iTrack2 = 0; iTrack2 < fVevent->GetNumberOfTracks(); iTrack2++)
 187     {
 188       if(iTrack1==iTrack2) continue;
 189
 190       AliVParticle* Vtrack2 = fVevent->GetTrack(iTrack2);
 191       if (!Vtrack2)
 192         {
 193           printf("ERROR: Could not receive track %d\n", iTrack2);
 194           continue;
 195         }
 196
 197       AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
 198       AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
 199       AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
 200
 201       if(fAlgorithm=="DCA")
 202       {
 203                   extTrackParam2 = new AliExternalTrackParam();
 204                         extTrackParam2->CopyFromVTrack(track2);
 205                 }
 206
 207       //Second track cuts
 208       if(fIsAOD)
 209         {
 210           if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
 211           if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
 212           if(atrack2->GetTPCNcls() < fTpcNcls) continue;
 213         }
 214       else
 215         {
 216           if(!fTrackCuts->AcceptTrack(etrack2)) continue;
 217         }
 218
 219       //Second track pid
 220       Double_t tpcNsigma2 = fPIDResponse->NumberOfSigmasTPC(track2,AliPID::kElectron);
 221       if(tpcNsigma2<fTPCnSigmaMin || tpcNsigma2>fTPCnSigmaMax) continue;
 222
 223       //Pt Cut
 224       if((track2->Pt() < fPtMin) && (fHasPtCut)) continue;
 225
 226       if(fAlgorithm=="DCA")
 227         {
 228           //Variables
 229           Double_t p1[3];
 230           Double_t p2[3];
 231           Double_t xt1; //radial position track 1 at the DCA point
 232           Double_t xt2; //radial position track 2 at the DCA point
 233           Double_t dca12; //DCA 1-2
 234           Bool_t hasdcaT1;
 235           Bool_t hasdcaT2;
 236
 237           if(fIsAOD)
 238           {
 239                   //DCA track1-track2
 240                   dca12 = extTrackParam2->GetDCA(extTrackParam1,bfield,xt2,xt1);
 241
 242                   //Momento of the track extrapolated to DCA track-track
 243                   //Track1
 244                   hasdcaT1 = extTrackParam1->GetPxPyPzAt(xt1,bfield,p1);
 245                   //Track2
 246                   hasdcaT2 = extTrackParam2->GetPxPyPzAt(xt2,bfield,p2);
 247           }
 248           else
 249           {
 250                   //DCA track1-track2
 251                   dca12 = etrack2->GetDCA(etrack1,bfield,xt2,xt1);
 252
 253                   //Momento of the track extrapolated to DCA track-track
 254                   //Track1
 255                   hasdcaT1 = etrack1->GetPxPyPzAt(xt1,bfield,p1);
 256                   //Track2
 257                   hasdcaT2 = etrack2->GetPxPyPzAt(xt2,bfield,p2);
 258           }
 259
 260           if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
 261
 262           //track1-track2 Invariant Mass
 263           Double_t eMass = 0.000510998910; //Electron mass in GeV
 264           Double_t pP1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]); //Track 1 momentum
 265           Double_t pP2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]); //Track 1 momentum
 266
 267           //Double_t E1 = sqrt(eMass*eMass+pP1*pP1);
 268           //Double_t E2 = sqrt(eMass*eMass+pP2*pP2);
 269           //Double_t imass = sqrt(2*eMass*eMass+2*(E1*E2-(p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2])));
 270           //Double_t angle = TMath::ACos((p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2])/(pP1*pP2));
 271
 272           TLorentzVector v1(p1[0],p1[1],p1[2],sqrt(eMass*eMass+pP1*pP1));
 273           TLorentzVector v2(p2[0],p2[1],p2[2],sqrt(eMass*eMass+pP2*pP2));
 274           Double_t imass = (v1+v2).M(); //Invariant Mass
 275           Double_t angle = v1.Angle(v2.Vect()); //Opening Angle (Total Angle)
 276
 277           Float_t fCharge1 = track1->Charge();
 278           Float_t fCharge2 = track2->Charge();
 279
 280           if(imass<fMassCut && angle<fAngleCut && dca12<fdcaCut)
 281             {
 282               if(fCharge1*fCharge2<0)
 283                 {
 284                   fIsULS=kTRUE;
 285                   fULSPartner[fNULS] = iTrack2;
 286                   fNULS++;
 287                 }
 288               if(fCharge1*fCharge2>0)
 289                 {
 290                   fIsLS=kTRUE;
 291                   fLSPartner[fNLS] = iTrack2;
 292                   fNLS++;
 293                 }
 294             }
 295
 296           //Fill some histograms
 297           if(fCharge1*fCharge2<0 && fHistMass) fHistMass->Fill(imass);
 298           if(fCharge1*fCharge2>0 && fHistMassBack) fHistMassBack->Fill(imass);
 299
 300           if(fCharge1*fCharge2<0 && fHistDCA) fHistDCA->Fill(dca12);
 301           if(fCharge1*fCharge2>0 && fHistDCABack) fHistDCABack->Fill(dca12);
 302
 303           if(fCharge1*fCharge2<0 && fHistAngle) fHistAngle->Fill(angle);
 304           if(fCharge1*fCharge2>0 && fHistAngleBack) fHistAngleBack->Fill(angle);
 305
 306         }
 307       else if(fAlgorithm=="KF")
 308         {
 309           Int_t fPDGtrack1 = 11;
 310           Int_t fPDGtrack2 = 11;
 311
 312           Float_t fCharge1 = track1->Charge();
 313           Float_t fCharge2 = track2->Charge();
 314
 315           if(fCharge1>0) fPDGtrack1 = -11;
 316           if(fCharge2>0) fPDGtrack2 = -11;
 317
 318           AliKFParticle::SetField(bfield);
 319           AliKFParticle fKFtrack1(*track1, fPDGtrack1);
 320           AliKFParticle fKFtrack2(*track2, fPDGtrack2);
 321           AliKFParticle fRecoGamma(fKFtrack1, fKFtrack2);
 322
 323           //Reconstruction Cuts
 324           if(fRecoGamma.GetNDF()<1) continue;
 325           Double_t chi2OverNDF = fRecoGamma.GetChi2()/fRecoGamma.GetNDF();
 326           if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
 327
 328           //Invariant Mass
 329           Double_t imass;
 330           Double_t width;
 331           fRecoGamma.GetMass(imass,width);
 332
 333           //Opening Angle (Total Angle)
 334           Double_t angle = fKFtrack1.GetAngle(fKFtrack2);
 335
 336           //Fill some histograms
 337           if(fCharge1*fCharge2<0 && fHistAngle) fHistAngle->Fill(angle);
 338           if(fCharge1*fCharge2>0 && fHistAngleBack) fHistAngleBack->Fill(angle);
 339
 340           if(angle>fAngleCut) continue;
 341
 342           if(imass<fMassCut)
 343             {
 344               if(fCharge1*fCharge2<0)
 345                 {
 346                   fIsULS=kTRUE;
 347                   fULSPartner[fNULS] = iTrack2;
 348                   fNULS++;
 349                 }
 350               if(fCharge1*fCharge2>0)
 351                 {
 352                   fIsLS=kTRUE;
 353                   fLSPartner[fNLS] = iTrack2;
 354                   fNLS++;
 355                 }
 356             }
 357
 358           //Fill some histograms
 359           if(fCharge1*fCharge2<0 && fHistMass) fHistMass->Fill(imass);
 360           if(fCharge1*fCharge2>0 && fHistMassBack) fHistMassBack->Fill(imass);
 361
 362         }
 363       else
 364         {
 365           AliError( Form("Error: %s is not a valid algorithm option.",(const char*)fAlgorithm));
 366           return;
 367         }
 368       if(extTrackParam2) delete extTrackParam2;
 369
 370     }
 371
 372
 373   if(extTrackParam1) delete extTrackParam1;
 374
 375
 376   return;
 377 }