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
192 // max nSig cuts for 10d
195 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
196 if(esdtrack==NULL)return maxCut;
197 double pt = esdtrack->Pt();
200 { //<--- new para for updated non-liniarity
201 double meanP[3] = {1.02902,1.22293,1.67287e-01};
202 double sigP[3] = {5.36355e-02,-3.25999e-02,4.42442e-01};
203 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
204 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
205 maxCut = mean+3.0*sig;
209 double meanP[3] = {0.99,1.299,0.35};
210 double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
211 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
212 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
213 maxCut = mean+3.0*sig;
218 Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
220 // mim nsig cuts for 10d
222 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
223 if(esdtrack==NULL)return mimCut;
224 double pt = esdtrack->Pt();
226 if(flageop<0.5) // real
229 //double meanP[3] = {1.056,1.169,0.1339};
230 //double sigP[3] = {3.37e-05,1.298e-04,1.403e-04};
231 double meanP[3] = {1.02902,1.22293,1.67287e-01};
232 double sigP[3] = {5.36355e-02,-3.25999e-02,4.42442e-01};
233 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
234 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
235 mimCut = mean+0.0*sig;
240 double meanP[3] = {0.99,1.299,0.35};
241 double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
242 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
243 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
244 mimCut = mean+0.0*sig;
249 //___________________________________________________________________________
250 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
251 { // Returns e/p & Rmatch
253 Double_t matchclsE = 9999.9;
254 double feop = -9999.9;
256 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
257 if(esdtrack==NULL)return feop;
258 const AliESDEvent *evt = esdtrack->GetESDEvent();
260 Int_t icl = esdtrack->GetEMCALcluster();
262 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
263 if(!cluster->IsEMCAL()) {return feop;}
265 matchclsE = cluster->E();
267 if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
275 //____________________________________________________________________________________
276 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
277 { // Returns e/p if an electron is matched
281 Double_t matchclsE = 9999.9;
283 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
284 AliESDEvent *evt = esdtrack->GetESDEvent();
285 Double_t magF = evt->GetMagneticField();
286 Double_t magSign = 1.0;
287 if(magF<0)magSign = -1.0;
288 //printf("magF ; %g ; %g \n", magF,magSign);
290 if (!TGeoGlobalMagField::Instance()->GetField()) {
291 printf("Loading field map...\n");
292 //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
293 AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
294 TGeoGlobalMagField::Instance()->SetField(field);
297 AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
299 emctrack->GetBxByBz(fieldB);
300 //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
302 for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
305 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
306 if(!cluster->IsEMCAL()) continue;
308 cluster->GetPosition(clsPos);
309 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
310 emctrack->GetXYZ(trkPos);
312 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
313 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
315 Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
316 Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
318 double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
322 matchclsE = cluster->E();
327 double feop = -9999.9;
328 if(matchclsE<9999) feop = matchclsE/esdtrack->P();
330 // if(feop!=-9999.9)printf("%f\n",feop) ;