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