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