Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliSelectNonHFE.cxx
CommitLineData
dfcc2025 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"
7edfba87 32#include "AliAODTrack.h"
dfcc2025 33#include "AliESDtrack.h"
7edfba87 34#include "AliVEvent.h"
35#include "AliESDtrackCuts.h"
36#include "AliVTrack.h"
37#include "AliVParticle.h"
70448e3b 38#include "AliPIDResponse.h"
39#include "AliPID.h"
7edfba87 40#include "AliExternalTrackParam.h"
dfcc2025 41#include "AliSelectNonHFE.h"
42#include "AliKFParticle.h"
43#include "AliLog.h"
44#include "stdio.h"
45#include "iostream"
46#include "fstream"
47
48ClassImp(AliSelectNonHFE)
49//________________________________________________________________________
50AliSelectNonHFE::AliSelectNonHFE(const char *name, const Char_t *title)
51: TNamed(name, title)
52 ,fTrackCuts(0)
7edfba87 53 ,fAlgorithm("DCA")
dfcc2025 54 ,fAngleCut(999)
55 ,fdcaCut(999)
70448e3b 56 ,fTPCnSigmaMin(-3)
57 ,fTPCnSigmaMax(3)
dfcc2025 58 ,fMassCut(0.5)
59 ,fChi2OverNDFCut(999)
a6f21076 60 ,fPtMin(0.3)
dfcc2025 61 ,fIsLS(kFALSE)
62 ,fIsULS(kFALSE)
7edfba87 63 ,fIsAOD(kFALSE)
a6f21076 64 ,fHasPtCut(kFALSE)
dfcc2025 65 ,fNLS(0)
66 ,fNULS(0)
a6f21076 67 ,fTpcNcls(50)
dfcc2025 68 ,fLSPartner(0)
69 ,fULSPartner(0)
70 ,fHistMass(0)
71 ,fHistMassBack(0)
72 ,fHistDCA(0)
73 ,fHistDCABack(0)
74 ,fHistAngle(0)
75 ,fHistAngleBack(0)
70448e3b 76 ,fPIDResponse(0)
dfcc2025 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//________________________________________________________________________
96AliSelectNonHFE::AliSelectNonHFE()
97 : TNamed()
98 ,fTrackCuts(0)
7edfba87 99 ,fAlgorithm("DCA")
dfcc2025 100 ,fAngleCut(999)
101 ,fdcaCut(999)
70448e3b 102 ,fTPCnSigmaMin(-3)
103 ,fTPCnSigmaMax(3)
dfcc2025 104 ,fMassCut(0.5)
105 ,fChi2OverNDFCut(999)
a6f21076 106 ,fPtMin(0.3)
dfcc2025 107 ,fIsLS(kFALSE)
108 ,fIsULS(kFALSE)
7edfba87 109 ,fIsAOD(kFALSE)
a6f21076 110 ,fHasPtCut(kFALSE)
dfcc2025 111 ,fNLS(0)
112 ,fNULS(0)
a6f21076 113 ,fTpcNcls(50)
dfcc2025 114 ,fLSPartner(0)
115 ,fULSPartner(0)
116 ,fHistMass(0)
117 ,fHistMassBack(0)
118 ,fHistDCA(0)
119 ,fHistDCABack(0)
120 ,fHistAngle(0)
121 ,fHistAngleBack(0)
70448e3b 122 ,fPIDResponse(0)
dfcc2025 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//_________________________________________
142AliSelectNonHFE::~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//__________________________________________
7edfba87 156void AliSelectNonHFE::FindNonHFE(Int_t iTrack1, AliVParticle *Vtrack1, AliVEvent *fVevent)
dfcc2025 157{
7edfba87 158 AliVTrack *track1 = dynamic_cast<AliVTrack*>(Vtrack1);
159 AliESDtrack *etrack1 = dynamic_cast<AliESDtrack*>(Vtrack1);
2d17adb0 160 AliExternalTrackParam *extTrackParam1=NULL;
161 AliExternalTrackParam *extTrackParam2=NULL;
7edfba87 162 if(fAlgorithm=="DCA")
163 {
164 extTrackParam1 = new AliExternalTrackParam();
165 extTrackParam1->CopyFromVTrack(track1);
166 }
dfcc2025 167 //
168 // Find non HFE electrons
169 //
170
dfcc2025 171 //Magnetic Field
7edfba87 172 Double_t bfield = fVevent->GetMagneticField();
dfcc2025 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
7edfba87 186 for(Int_t iTrack2 = 0; iTrack2 < fVevent->GetNumberOfTracks(); iTrack2++)
dfcc2025 187 {
7edfba87 188 if(iTrack1==iTrack2) continue;
7c4ec6e7 189
7edfba87 190 AliVParticle* Vtrack2 = fVevent->GetTrack(iTrack2);
191 if (!Vtrack2)
dfcc2025 192 {
193 printf("ERROR: Could not receive track %d\n", iTrack2);
194 continue;
195 }
7edfba87 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
dfcc2025 207 //Second track cuts
7edfba87 208 if(fIsAOD)
209 {
210 if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
211 if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
a6f21076 212 if(atrack2->GetTPCNcls() < fTpcNcls) continue;
7edfba87 213 }
214 else
215 {
216 if(!fTrackCuts->AcceptTrack(etrack2)) continue;
217 }
218
219 //Second track pid
70448e3b 220 Double_t tpcNsigma2 = fPIDResponse->NumberOfSigmasTPC(track2,AliPID::kElectron);
70448e3b 221 if(tpcNsigma2<fTPCnSigmaMin || tpcNsigma2>fTPCnSigmaMax) continue;
dfcc2025 222
a6f21076 223 //Pt Cut
224 if((track2->Pt() < fPtMin) && (fHasPtCut)) continue;
7edfba87 225
226 if(fAlgorithm=="DCA")
dfcc2025 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
7edfba87 233 Double_t dca12; //DCA 1-2
234 Bool_t hasdcaT1;
235 Bool_t hasdcaT2;
dfcc2025 236
7edfba87 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 }
dfcc2025 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
a6f21076 318 AliKFParticle::SetField(bfield);
dfcc2025 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
7c4ec6e7 340 if(angle>fAngleCut) continue;
dfcc2025 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 }
2d17adb0 368 if(extTrackParam2) delete extTrackParam2;
369
dfcc2025 370 }
7edfba87 371
372
2d17adb0 373 if(extTrackParam1) delete extTrackParam1;
374
dfcc2025 375
376 return;
377}