1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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)
27 #include "AliESDInputHandler.h"
28 #include "AliESDpid.h"
30 #include "AliHFEdetPIDqa.h"
31 #include "AliHFEpidEMCAL.h"
32 #include "AliHFEpidQAmanager.h"
34 #include "AliHFEemcalPIDqa.h"
35 //#include "AliVCluster.h"
36 //#include "AliVCaloCells.h"
37 //#include "AliVEvent.h"
40 //#include "AliESDEvent.h"
41 //#include "AliESDtrack.h"
42 #include "AliHFEpidEMCAL.h"
44 ClassImp(AliHFEpidEMCAL)
46 //___________________________________________________________________
47 AliHFEpidEMCAL::AliHFEpidEMCAL():
58 //___________________________________________________________________
59 AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
70 //___________________________________________________________________
71 AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
83 //___________________________________________________________________
84 AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
86 // Assignment operator
95 //___________________________________________________________________
96 AliHFEpidEMCAL::~AliHFEpidEMCAL(){
100 if(fPID) delete fPID;
102 //___________________________________________________________________
103 void AliHFEpidEMCAL::Copy(TObject &ref) const {
105 // Performs the copying of the object
107 AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
110 target.feopMax = feopMax;
111 target.feopMim = feopMim;
113 AliHFEpidBase::Copy(ref);
115 //___________________________________________________________________
116 Bool_t AliHFEpidEMCAL::InitializePID(Int_t /*run*/){
118 // InitializePID: EMCAL experts have to implement code here
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
128 if(track==NULL)return 0;
130 if(!fkPIDResponse) return 0;
131 AliDebug(2, "PID object available");
132 // EMCal not fESDpid (s.s Feb. 11)
134 const AliVTrack *trk = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
135 if (trk == NULL) return 0;
136 //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
137 if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
138 AliDebug(2, "Track Has EMCAL PID");
141 pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
142 // not QA for EMCal, will be added (s.s Feb. 11)
145 double feopMimCut = 0.0;
146 double feopMaxCut = 0.0;
148 feopMimCut = feopMim;
149 feopMaxCut = feopMax;
151 if(feopMax >900) // not use fix cut
153 feopMimCut = CalEopCutMim(track->GetRecTrack(),0);
154 feopMaxCut = CalEopCutMax(track->GetRecTrack(),0);
156 if(feopMax < -900) // not use fix cut for MC
158 feopMimCut = CalEopCutMim(track->GetRecTrack(),1);
159 feopMaxCut = CalEopCutMax(track->GetRecTrack(),1);
162 //printf("eop cuts org; %g; %g \n",feopMim,feopMax);
163 //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
165 Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
166 AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
167 AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax));
169 //if(eop>feopMim && eop<feopMax){
170 //if(eop>feopMim && eop<feopMax){
171 if(eop>feopMimCut && eop<feopMaxCut){
173 //printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop);
175 pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
179 pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
182 AliDebug(1, Form("eID %g ; %d \n",eop,pdg));
183 //printf("eID %g ; %d \n",eop,pdg);
189 //__________________________________________________________________________
190 Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const
193 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
194 if(esdtrack==NULL)return maxCut;
195 double Pt = esdtrack->Pt();
199 double meanP[3] = {0.991,1.0819,0.235};
200 double sigP[3] = {6.52e-02,2.04e-01,3.34e-01};
201 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt);
202 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
203 maxCut = mean+3.0*sig;
207 double meanP[3] = {0.99,1.299,0.35};
208 double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
209 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt);
210 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
211 maxCut = mean+3.0*sig;
216 Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
219 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
220 if(esdtrack==NULL)return mimCut;
221 double Pt = esdtrack->Pt();
223 if(flageop<0.5) // real
226 double meanP[3] = {0.991,1.0819,0.235};
227 double sigP[3] = {6.52e-02,2.04e-01,3.34e-01};
228 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt);
229 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
230 mimCut = mean-2.0*sig;
235 double meanP[3] = {0.99,1.299,0.35};
236 double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
237 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt);
238 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
239 mimCut = mean-2.0*sig;
244 //___________________________________________________________________________
245 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
246 { // Returns e/p & Rmatch
248 Double_t matchclsE = 9999.9;
249 double feop = -9999.9;
251 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
252 if(esdtrack==NULL)return feop;
253 const AliESDEvent *evt = esdtrack->GetESDEvent();
255 Int_t icl = esdtrack->GetEMCALcluster();
257 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
258 if(!cluster->IsEMCAL()) {return feop;}
260 matchclsE = cluster->E();
263 if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
271 //____________________________________________________________________________________
272 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
273 { // Returns e/p if an electron is matched
277 Double_t matchclsE = 9999.9;
279 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
280 AliESDEvent *evt = esdtrack->GetESDEvent();
281 Double_t magF = evt->GetMagneticField();
282 Double_t magSign = 1.0;
283 if(magF<0)magSign = -1.0;
284 //printf("magF ; %g ; %g \n", magF,magSign);
286 if (!TGeoGlobalMagField::Instance()->GetField()) {
287 printf("Loading field map...\n");
288 //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
289 AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
290 TGeoGlobalMagField::Instance()->SetField(field);
293 AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
295 emctrack->GetBxByBz(fieldB);
296 //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
298 for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
301 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
302 if(!cluster->IsEMCAL()) continue;
304 cluster->GetPosition(clsPos);
305 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
306 emctrack->GetXYZ(trkPos);
308 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
309 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
311 Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
312 Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
314 double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
318 matchclsE = cluster->E();
323 double feop = -9999.9;
324 if(matchclsE<9999) feop = matchclsE/esdtrack->P();
326 // if(feop!=-9999.9)printf("%f\n",feop) ;