]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGHF/hfe/AliSelectNonHFE.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / 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 }