Deepa code electron hadron correlation
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidEMCAL.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 // Class for EMCAL PID
17 // Implements the abstract base class AliHFEpidBase
18 // IsInitialized() does the PID decision with energy and 
19 // momentum matching (e/p) 
20 // 
21 // Authors:
22 //  Shingo Sakai 
23 //   
24 //
25 //
26 #include <TMath.h>
27 #include <TRandom3.h>
28 #include "AliESDInputHandler.h"
29 #include "AliESDpid.h"
30
31 #include "AliHFEdetPIDqa.h"
32 #include "AliHFEpidEMCAL.h"
33 #include "AliHFEpidQAmanager.h"
34
35 #include "AliHFEemcalPIDqa.h"
36 //#include "AliVCluster.h"
37 //#include "AliVCaloCells.h"
38 //#include "AliVEvent.h"
39 #include "AliLog.h"
40 #include "AliPID.h"
41 //#include "AliESDEvent.h"
42 //#include "AliESDtrack.h"
43 #include "AliHFEpidEMCAL.h"
44
45 ClassImp(AliHFEpidEMCAL)
46
47 //___________________________________________________________________
48 AliHFEpidEMCAL::AliHFEpidEMCAL():
49   AliHFEpidBase()
50   , fPID(NULL)
51   , feopMim(0.8)
52   , feopMax(1.4)
53 {
54   //
55   // Constructor
56   //
57
58
59 //___________________________________________________________________
60 AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
61   AliHFEpidBase(name)
62   , fPID(NULL)
63   , feopMim(0.8)
64   , feopMax(1.4)
65 {
66   //
67   // Constructor
68   //
69 }
70
71 //___________________________________________________________________
72 AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
73   AliHFEpidBase("")
74   , fPID(NULL)
75   , feopMim(0.8)
76   , feopMax(1.4)
77 {  
78   // 
79   // Copy operator
80   //
81
82   c.Copy(*this);
83 }
84 //___________________________________________________________________
85 AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
86   //
87   // Assignment operator
88   //
89
90   if(this != &ref){
91     ref.Copy(*this);
92   }
93
94   return *this;
95 }
96 //___________________________________________________________________
97 AliHFEpidEMCAL::~AliHFEpidEMCAL(){
98   //
99   // Destructor
100   //
101   if(fPID) delete fPID;
102 }
103 //___________________________________________________________________
104 void AliHFEpidEMCAL::Copy(TObject &ref) const {
105   //
106   // Performs the copying of the object
107   //
108   AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
109
110   target.fPID = fPID;          
111   target.feopMax = feopMax;
112   target.feopMim = feopMim;
113
114   AliHFEpidBase::Copy(ref);
115 }
116 //___________________________________________________________________
117 Bool_t AliHFEpidEMCAL::InitializePID(Int_t /*run*/){
118   //
119   // InitializePID: EMCAL experts have to implement code here
120   //
121   return kTRUE;
122 }
123
124 //___________________________________________________________________
125 Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
126 { //Function to a return a code indicating whether or not an electron candidate is selected
127   //
128   //
129   if(track==NULL)return 0;
130
131   if(!fkPIDResponse) return 0;
132   AliDebug(2, "PID object available");
133   // EMCal not fESDpid  (s.s Feb. 11)
134         
135   const AliVTrack *trk = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
136   if (trk == NULL) return 0;
137   //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
138   if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
139   AliDebug(2, "Track Has EMCAL PID");
140   
141   if(pidqa) 
142   pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
143   // not QA for EMCal, will be added (s.s Feb. 11)
144
145
146      double feopMimCut = 0.0; 
147      double feopMaxCut = 0.0; 
148
149      feopMimCut = feopMim;
150      feopMaxCut = feopMax;
151
152      if(feopMax >900)  // not use fix cut
153        {
154         feopMimCut = CalEopCutMim(track->GetRecTrack(),0);
155         feopMaxCut = CalEopCutMax(track->GetRecTrack(),0); 
156        }
157      if(feopMax < -900)  // not use fix cut for MC
158        {
159         feopMimCut = CalEopCutMim(track->GetRecTrack(),1);
160         feopMaxCut = CalEopCutMax(track->GetRecTrack(),1); 
161        }
162
163   //printf("eop cuts org; %g; %g \n",feopMim,feopMax);
164   //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
165
166   Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
167   AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
168   AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax));
169   Int_t pdg = 0;
170   //if(eop>feopMim && eop<feopMax){
171   //if(eop>feopMim && eop<feopMax){
172   if(eop>feopMimCut && eop<feopMaxCut){
173     pdg = 11;
174     //printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop);
175     if(pidqa) 
176     pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
177   }
178   else
179   {
180    pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
181   } 
182
183     AliDebug(1, Form("eID %g ; %d \n",eop,pdg));  
184     //printf("eID %g ; %d \n",eop,pdg);  
185
186   return pdg;
187 }
188
189
190 //__________________________________________________________________________
191 Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const
192 {
193  // max nSig cuts for 10d
194
195   double maxCut = 0.0;
196   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
197   if(esdtrack==NULL)return maxCut;
198   double pt = esdtrack->Pt();
199
200   if(flageop<0.5)
201     { //<--- new para for updated non-liniarity
202      double meanP[3] = {1.056,1.169,0.1339}; 
203      double sigP[3] = {3.37e-05,1.298e-04,1.403e-04}; 
204      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
205      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
206      maxCut = mean+3.0*sig; 
207     }
208   else
209     {
210      double meanP[3] = {0.99,1.299,0.35}; 
211      double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; 
212      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
213      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
214      maxCut = mean+3.0*sig; 
215     }
216      return maxCut;
217 }
218
219 Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
220 {
221  // mim nsig cuts for 10d
222   double mimCut = 0.0;
223   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
224   if(esdtrack==NULL)return mimCut;
225   double pt = esdtrack->Pt();
226
227   if(flageop<0.5) // real
228      { 
229       //printf("real"); 
230      double meanP[3] = {1.056,1.169,0.1339}; 
231      double sigP[3] = {3.37e-05,1.298e-04,1.403e-04}; 
232      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
233      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
234      mimCut = mean-2.0*sig;
235     }
236    else // MC
237     { 
238       //printf("MC"); 
239      double meanP[3] = {0.99,1.299,0.35}; 
240      double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; 
241      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
242      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
243      mimCut = mean-2.0*sig;
244     }
245      return mimCut;
246 }
247
248 //___________________________________________________________________________
249 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
250 { // Returns e/p & Rmatch
251
252   Double_t matchclsE = 9999.9;
253   double feop = -9999.9;
254
255   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
256   if(esdtrack==NULL)return feop;
257   const AliESDEvent *evt = esdtrack->GetESDEvent();
258
259    Int_t icl = esdtrack->GetEMCALcluster();
260
261    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
262    if(!cluster->IsEMCAL()) {return feop;}
263
264        matchclsE = cluster->E();
265   
266    //printf("energy org  %f \n",matchclsE);
267    //printf("Label %d \n",esdtrack->GetLabel());
268    if(esdtrack->GetLabel()>-1)
269      {
270             double energyMC = cluster->E();
271             TRandom3 *smear = new TRandom3(icl); 
272             matchclsE = smear->Gaus(energyMC,0.07 * TMath::Sqrt(energyMC));
273      } 
274
275    //printf("energy af  %f \n",matchclsE);
276
277    if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
278
279    return feop;
280
281 }
282
283
284 /*
285 //____________________________________________________________________________________
286 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
287 { // Returns e/p if an electron is matched
288
289   Float_t  clsPos[3];
290   Double_t trkPos[3];
291   Double_t matchclsE = 9999.9;
292
293   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
294   AliESDEvent *evt = esdtrack->GetESDEvent();
295   Double_t  magF = evt->GetMagneticField();
296   Double_t magSign = 1.0;
297   if(magF<0)magSign = -1.0;
298   //printf("magF ; %g ; %g \n", magF,magSign);
299  
300   if (!TGeoGlobalMagField::Instance()->GetField()) {
301           printf("Loading field map...\n");
302           //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
303           AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
304           TGeoGlobalMagField::Instance()->SetField(field);
305                     }
306
307   AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
308   Double_t fieldB[3]; 
309   emctrack->GetBxByBz(fieldB);
310   //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
311
312   for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
313    {
314
315    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
316    if(!cluster->IsEMCAL()) continue;
317
318    cluster->GetPosition(clsPos);
319    if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
320    emctrack->GetXYZ(trkPos);
321
322    TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
323    TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
324
325    Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
326    Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
327
328    double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
329
330    if(rmatch<0.02)
331       {
332        matchclsE = cluster->E();
333       }
334   }
335    delete emctrack;
336
337    double feop = -9999.9;
338    if(matchclsE<9999) feop = matchclsE/esdtrack->P();
339
340    //   if(feop!=-9999.9)printf("%f\n",feop) ; 
341    return feop;
342
343 }
344 */
345