Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / 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 "AliESDInputHandler.h"
28 #include "AliESDpid.h"
29
30 #include "AliHFEdetPIDqa.h"
31 #include "AliHFEpidEMCAL.h"
32 #include "AliHFEpidQAmanager.h"
33
34 #include "AliHFEemcalPIDqa.h"
35 //#include "AliVCluster.h"
36 //#include "AliVCaloCells.h"
37 //#include "AliVEvent.h"
38 #include "AliLog.h"
39 #include "AliPID.h"
40 //#include "AliESDEvent.h"
41 //#include "AliESDtrack.h"
42 #include "AliHFEpidEMCAL.h"
43
44 ClassImp(AliHFEpidEMCAL)
45
46 //___________________________________________________________________
47 AliHFEpidEMCAL::AliHFEpidEMCAL():
48   AliHFEpidBase()
49   , fPID(NULL)
50   , feopMim(0.9)
51   , feopMax(1.3)
52 {
53   //
54   // Constructor
55   //
56
57
58 //___________________________________________________________________
59 AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
60   AliHFEpidBase(name)
61   , fPID(NULL)
62   , feopMim(0.9)
63   , feopMax(1.3)
64 {
65   //
66   // Constructor
67   //
68 }
69
70 //___________________________________________________________________
71 AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
72   AliHFEpidBase("")
73   , fPID(NULL)
74   , feopMim(0.9)
75   , feopMax(1.3)
76 {  
77   // 
78   // Copy operator
79   //
80
81   c.Copy(*this);
82 }
83 //___________________________________________________________________
84 AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
85   //
86   // Assignment operator
87   //
88
89   if(this != &ref){
90     ref.Copy(*this);
91   }
92
93   return *this;
94 }
95 //___________________________________________________________________
96 AliHFEpidEMCAL::~AliHFEpidEMCAL(){
97   //
98   // Destructor
99   //
100   if(fPID) delete fPID;
101 }
102 //___________________________________________________________________
103 void AliHFEpidEMCAL::Copy(TObject &ref) const {
104   //
105   // Performs the copying of the object
106   //
107   AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
108
109   target.fPID = fPID;          
110   target.feopMax = feopMax;
111   target.feopMim = feopMim;
112
113   AliHFEpidBase::Copy(ref);
114 }
115 //___________________________________________________________________
116 Bool_t AliHFEpidEMCAL::InitializePID(Int_t /*run*/){
117   //
118   // InitializePID: EMCAL experts have to implement code here
119   //
120   return kTRUE;
121 }
122
123 //___________________________________________________________________
124 Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
125 { //Function to a return a code indicating whether or not an electron candidate is selected
126   //
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.08833e+00,1.17837e+00,7.47923e-02}; 
203      double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-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] = {9.94183e-01,1.23952e+00,3.55571e-01}; 
211      double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-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
217   //printf("eop cuts org; %g; %g ; %g\n",feopMim,feopMax,pt);
218
219      return maxCut;
220 }
221
222 Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
223 {
224  // mim nsig cuts for 10d
225   double mimCut = 0.0;
226   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
227   if(esdtrack==NULL)return mimCut;
228   double pt = esdtrack->Pt();
229
230   if(flageop<0.5) // real
231      { 
232       //printf("real"); 
233      double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02}; 
234      double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; 
235      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
236      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
237      mimCut = mean-3.0*sig;
238     }
239    else // MC
240     { 
241       //printf("MC"); 
242      double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01}; 
243      double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01}; 
244      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
245      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
246      mimCut = mean-3.0*sig;
247     }
248      return mimCut;
249 }
250
251 //___________________________________________________________________________
252 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
253 { // Returns e/p & Rmatch
254
255   Double_t matchclsE = 9999.9;
256   double feop = -9999.9;
257
258   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
259   if(esdtrack==NULL)return feop;
260   const AliESDEvent *evt = esdtrack->GetESDEvent();
261
262    Int_t icl = esdtrack->GetEMCALcluster();
263
264    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
265    if(cluster==NULL) {return feop;}
266    if(!cluster->IsEMCAL()) {return feop;}
267
268        matchclsE = cluster->E();
269
270    if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
271
272    return feop;
273
274 }
275
276
277 /*
278 //____________________________________________________________________________________
279 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
280 { // Returns e/p if an electron is matched
281
282   Float_t  clsPos[3];
283   Double_t trkPos[3];
284   Double_t matchclsE = 9999.9;
285
286   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
287   AliESDEvent *evt = esdtrack->GetESDEvent();
288   Double_t  magF = evt->GetMagneticField();
289   Double_t magSign = 1.0;
290   if(magF<0)magSign = -1.0;
291   //printf("magF ; %g ; %g \n", magF,magSign);
292  
293   if (!TGeoGlobalMagField::Instance()->GetField()) {
294           printf("Loading field map...\n");
295           //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
296           AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
297           TGeoGlobalMagField::Instance()->SetField(field);
298                     }
299
300   AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
301   Double_t fieldB[3]; 
302   emctrack->GetBxByBz(fieldB);
303   //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
304
305   for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
306    {
307
308    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
309    if(!cluster->IsEMCAL()) continue;
310
311    cluster->GetPosition(clsPos);
312    if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
313    emctrack->GetXYZ(trkPos);
314
315    TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
316    TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
317
318    Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
319    Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
320
321    double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
322
323    if(rmatch<0.02)
324       {
325        matchclsE = cluster->E();
326       }
327   }
328    delete emctrack;
329
330    double feop = -9999.9;
331    if(matchclsE<9999) feop = matchclsE/esdtrack->P();
332
333    //   if(feop!=-9999.9)printf("%f\n",feop) ; 
334    return feop;
335
336 }
337 */
338