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
  33 #include "AliESDEvent.h"
  34 #include "AliESDtrackCuts.h"
  35 #include "AliESDtrack.h"
  36
  37 #include "AliSelectNonHFE.h"
  38 #include "AliKFParticle.h"
  39 #include "AliLog.h"
  40 #include "stdio.h"
  41 #include "iostream"
  42 #include "fstream"
  43
  44 ClassImp(AliSelectNonHFE)
  45 //________________________________________________________________________
  46 AliSelectNonHFE::AliSelectNonHFE(const char *name, const Char_t *title)
  47 : TNamed(name, title)
  48   ,fTrackCuts(0)
  49   ,fAlgorithm("MA")
  50   ,fAngleCut(999)
  51   ,fdcaCut(999)
  52   ,fdEdxMin(62)
  53   ,fdEdxMax(100)
  54   ,fMassCut(0.5)
  55   ,fChi2OverNDFCut(999)
  56   ,fIsLS(kFALSE)
  57   ,fIsULS(kFALSE)
  58   ,fNLS(0)
  59   ,fNULS(0)
  60   ,fLSPartner(0)
  61   ,fULSPartner(0)
  62   ,fHistMass(0)
  63   ,fHistMassBack(0)
  64   ,fHistDCA(0)
  65   ,fHistDCABack(0)
  66   ,fHistAngle(0)
  67   ,fHistAngleBack(0)
  68 {
  69   //
  70   // Constructor
  71   //
  72
  73   fTrackCuts = new AliESDtrackCuts();
  74   //Configure Default Track Cuts
  75   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  76   fTrackCuts->SetRequireTPCRefit(kTRUE);
  77   fTrackCuts->SetEtaRange(-0.9,0.9);
  78   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
  79   fTrackCuts->SetMaxChi2PerClusterTPC(4.0);
  80   fTrackCuts->SetMinNClustersTPC(50);
  81   fTrackCuts->SetPtRange(0.3,1e10);
  82
  83
  84 }
  85
  86 //________________________________________________________________________
  87 AliSelectNonHFE::AliSelectNonHFE()
  88   : TNamed()
  89   ,fTrackCuts(0)
  90   ,fAlgorithm("MA")
  91   ,fAngleCut(999)
  92   ,fdcaCut(999)
  93   ,fdEdxMin(62)
  94   ,fdEdxMax(100)
  95   ,fMassCut(0.5)
  96   ,fChi2OverNDFCut(999)
  97   ,fIsLS(kFALSE)
  98   ,fIsULS(kFALSE)
  99   ,fNLS(0)
 100   ,fNULS(0)
 101   ,fLSPartner(0)
 102   ,fULSPartner(0)
 103   ,fHistMass(0)
 104   ,fHistMassBack(0)
 105   ,fHistDCA(0)
 106   ,fHistDCABack(0)
 107   ,fHistAngle(0)
 108   ,fHistAngleBack(0)
 109 {
 110   //
 111   // Constructor
 112   //
 113
 114   fTrackCuts = new AliESDtrackCuts();
 115   //Configure Default Track Cuts
 116   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 117   fTrackCuts->SetRequireTPCRefit(kTRUE);
 118   fTrackCuts->SetEtaRange(-0.9,0.9);
 119   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 120   fTrackCuts->SetMaxChi2PerClusterTPC(4.0);
 121   fTrackCuts->SetMinNClustersTPC(50);
 122   fTrackCuts->SetPtRange(0.3,1e10);
 123
 124
 125 }
 126
 127 //_________________________________________
 128 AliSelectNonHFE::~AliSelectNonHFE()
 129 {
 130   //
 131   // Destructor
 132   //
 133
 134   if(fTrackCuts) delete fTrackCuts;
 135   if(fLSPartner) delete [] fLSPartner;
 136   if(fULSPartner) delete [] fULSPartner;
 137
 138
 139 }
 140
 141 //__________________________________________
 142 void AliSelectNonHFE::FindNonHFE(Int_t iTrack1, AliESDtrack *track1, AliESDEvent *fESD)
 143 {
 144   //
 145   // Find non HFE electrons
 146   //
 147
 148
 149   //Magnetic Field
 150   Double_t bfield = fESD->GetMagneticField();
 151
 152   //Second Track loop
 153
 154   fIsULS = kFALSE;      //Non-HFE Unlike signal Flag
 155   fIsLS = kFALSE;       //Non-HFE like signal Flag
 156   fNULS = 0;            //Non-HFE Unlike signal Flag
 157   fNLS = 0;             //Non-HFE like signal Flag
 158
 159   if(fLSPartner) delete [] fLSPartner;
 160   if(fULSPartner) delete [] fULSPartner;
 161   fLSPartner = new int [100];   //store the partners index
 162   fULSPartner = new int [100];  //store the partners index
 163
 164   for(Int_t iTrack2 = iTrack1+1; iTrack2 < fESD->GetNumberOfTracks(); iTrack2++)
 165     {
 166       AliESDtrack* track2 = fESD->GetTrack(iTrack2);
 167       if (!track2)
 168         {
 169           printf("ERROR: Could not receive track %d\n", iTrack2);
 170           continue;
 171         }
 172
 173       //Second track cuts
 174       Double_t dEdx2 = track2->GetTPCsignal();
 175       if(dEdx2<fdEdxMin || dEdx2>fdEdxMax) continue;
 176       if(!fTrackCuts->AcceptTrack(track2)) continue;
 177
 178       if(fAlgorithm=="MA")
 179         {
 180           //Variables
 181           Double_t p1[3];
 182           Double_t p2[3];
 183           Double_t xt1; //radial position track 1 at the DCA point
 184           Double_t xt2; //radial position track 2 at the DCA point
 185           //DCA track1-track2
 186           Double_t dca12 = track2->GetDCA(track1,bfield,xt2,xt1);
 187
 188           //Momento of the track extrapolated to DCA track-track
 189           //Track1
 190           Bool_t hasdcaT1 = track1->GetPxPyPzAt(xt1,bfield,p1);
 191           //Track2
 192           Bool_t hasdcaT2 = track2->GetPxPyPzAt(xt2,bfield,p2);
 193
 194           if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
 195
 196           //track1-track2 Invariant Mass
 197           Double_t eMass = 0.000510998910; //Electron mass in GeV
 198           Double_t pP1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]); //Track 1 momentum
 199           Double_t pP2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]); //Track 1 momentum
 200
 201           //Double_t E1 = sqrt(eMass*eMass+pP1*pP1);
 202           //Double_t E2 = sqrt(eMass*eMass+pP2*pP2);
 203           //Double_t imass = sqrt(2*eMass*eMass+2*(E1*E2-(p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2])));
 204           //Double_t angle = TMath::ACos((p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2])/(pP1*pP2));
 205
 206           TLorentzVector v1(p1[0],p1[1],p1[2],sqrt(eMass*eMass+pP1*pP1));
 207           TLorentzVector v2(p2[0],p2[1],p2[2],sqrt(eMass*eMass+pP2*pP2));
 208           Double_t imass = (v1+v2).M(); //Invariant Mass
 209           Double_t angle = v1.Angle(v2.Vect()); //Opening Angle (Total Angle)
 210
 211           Float_t fCharge1 = track1->Charge();
 212           Float_t fCharge2 = track2->Charge();
 213
 214           if(imass<fMassCut && angle<fAngleCut && dca12<fdcaCut)
 215             {
 216               if(fCharge1*fCharge2<0)
 217                 {
 218                   fIsULS=kTRUE;
 219                   fULSPartner[fNULS] = iTrack2;
 220                   fNULS++;
 221                 }
 222               if(fCharge1*fCharge2>0)
 223                 {
 224                   fIsLS=kTRUE;
 225                   fLSPartner[fNLS] = iTrack2;
 226                   fNLS++;
 227                 }
 228             }
 229
 230           //Fill some histograms
 231           if(fCharge1*fCharge2<0 && fHistMass) fHistMass->Fill(imass);
 232           if(fCharge1*fCharge2>0 && fHistMassBack) fHistMassBack->Fill(imass);
 233
 234           if(fCharge1*fCharge2<0 && fHistDCA) fHistDCA->Fill(dca12);
 235           if(fCharge1*fCharge2>0 && fHistDCABack) fHistDCABack->Fill(dca12);
 236
 237           if(fCharge1*fCharge2<0 && fHistAngle) fHistAngle->Fill(angle);
 238           if(fCharge1*fCharge2>0 && fHistAngleBack) fHistAngleBack->Fill(angle);
 239
 240         }
 241       else if(fAlgorithm=="KF")
 242         {
 243           Int_t fPDGtrack1 = 11;
 244           Int_t fPDGtrack2 = 11;
 245
 246           Float_t fCharge1 = track1->Charge();
 247           Float_t fCharge2 = track2->Charge();
 248
 249           if(fCharge1>0) fPDGtrack1 = -11;
 250           if(fCharge2>0) fPDGtrack2 = -11;
 251
 252           AliKFParticle fKFtrack1(*track1, fPDGtrack1);
 253           AliKFParticle fKFtrack2(*track2, fPDGtrack2);
 254           AliKFParticle fRecoGamma(fKFtrack1, fKFtrack2);
 255
 256           //Reconstruction Cuts
 257           if(fRecoGamma.GetNDF()<1) continue;
 258           Double_t chi2OverNDF = fRecoGamma.GetChi2()/fRecoGamma.GetNDF();
 259           if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
 260
 261           //Invariant Mass
 262           Double_t imass;
 263           Double_t width;
 264           fRecoGamma.GetMass(imass,width);
 265
 266           //Opening Angle (Total Angle)
 267           Double_t angle = fKFtrack1.GetAngle(fKFtrack2);
 268
 269           //Fill some histograms
 270           if(fCharge1*fCharge2<0 && fHistAngle) fHistAngle->Fill(angle);
 271           if(fCharge1*fCharge2>0 && fHistAngleBack) fHistAngleBack->Fill(angle);
 272
 273           if(angle<fAngleCut) continue;
 274
 275           if(imass<fMassCut)
 276             {
 277               if(fCharge1*fCharge2<0)
 278                 {
 279                   fIsULS=kTRUE;
 280                   fULSPartner[fNULS] = iTrack2;
 281                   fNULS++;
 282                 }
 283               if(fCharge1*fCharge2>0)
 284                 {
 285                   fIsLS=kTRUE;
 286                   fLSPartner[fNLS] = iTrack2;
 287                   fNLS++;
 288                 }
 289             }
 290
 291           //Fill some histograms
 292           if(fCharge1*fCharge2<0 && fHistMass) fHistMass->Fill(imass);
 293           if(fCharge1*fCharge2>0 && fHistMassBack) fHistMassBack->Fill(imass);
 294
 295         }
 296       else
 297         {
 298           AliError( Form("Error: %s is not a valid algorithm option.",(const char*)fAlgorithm));
 299           return;
 300         }
 301
 302     }
 303
 304   return;
 305 }