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.056,1.169,0.1339};
202 double sigP[3] = {3.37e-05,1.298e-04,1.403e-04};
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 mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
232 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
233 mimCut = mean-2.0*sig;
238 double meanP[3] = {0.99,1.299,0.35};
239 double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
240 double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
241 double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
242 mimCut = mean-2.0*sig;
247 //___________________________________________________________________________
248 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
249 { // Returns e/p & Rmatch
251 Double_t matchclsE = 9999.9;
252 double feop = -9999.9;
254 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
255 if(esdtrack==NULL)return feop;
256 const AliESDEvent *evt = esdtrack->GetESDEvent();
258 Int_t icl = esdtrack->GetEMCALcluster();
260 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
261 if(!cluster->IsEMCAL()) {return feop;}
263 matchclsE = cluster->E();
265 if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
273 //____________________________________________________________________________________
274 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
275 { // Returns e/p if an electron is matched
279 Double_t matchclsE = 9999.9;
281 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
282 AliESDEvent *evt = esdtrack->GetESDEvent();
283 Double_t magF = evt->GetMagneticField();
284 Double_t magSign = 1.0;
285 if(magF<0)magSign = -1.0;
286 //printf("magF ; %g ; %g \n", magF,magSign);
288 if (!TGeoGlobalMagField::Instance()->GetField()) {
289 printf("Loading field map...\n");
290 //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
291 AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
292 TGeoGlobalMagField::Instance()->SetField(field);
295 AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
297 emctrack->GetBxByBz(fieldB);
298 //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
300 for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
303 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
304 if(!cluster->IsEMCAL()) continue;
306 cluster->GetPosition(clsPos);
307 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
308 emctrack->GetXYZ(trkPos);
310 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
311 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
313 Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
314 Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
316 double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
320 matchclsE = cluster->E();
325 double feop = -9999.9;
326 if(matchclsE<9999) feop = matchclsE/esdtrack->P();
328 // if(feop!=-9999.9)printf("%f\n",feop) ;