]>
Commit | Line | Data |
---|---|---|
c2690925 | 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" | |
c2690925 | 28 | #include "AliESDpid.h" |
29 | ||
30 | #include "AliHFEdetPIDqa.h" | |
e156c3bb | 31 | #include "AliHFEpidEMCAL.h" |
c2690925 | 32 | #include "AliHFEpidQAmanager.h" |
33 | ||
e156c3bb | 34 | #include "AliHFEemcalPIDqa.h" |
35 | //#include "AliVCluster.h" | |
36 | //#include "AliVCaloCells.h" | |
37 | //#include "AliVEvent.h" | |
c2690925 | 38 | #include "AliLog.h" |
c2690925 | 39 | #include "AliPID.h" |
e156c3bb | 40 | //#include "AliESDEvent.h" |
41 | //#include "AliESDtrack.h" | |
c2690925 | 42 | #include "AliHFEpidEMCAL.h" |
43 | ||
44 | ClassImp(AliHFEpidEMCAL) | |
45 | ||
46 | //___________________________________________________________________ | |
47 | AliHFEpidEMCAL::AliHFEpidEMCAL(): | |
48 | AliHFEpidBase() | |
49 | , fPID(NULL) | |
fb87d707 | 50 | , feopMim(0.9) |
51 | , feopMax(1.3) | |
c2690925 | 52 | { |
53 | // | |
54 | // Constructor | |
55 | // | |
56 | } | |
57 | ||
58 | //___________________________________________________________________ | |
59 | AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name): | |
60 | AliHFEpidBase(name) | |
61 | , fPID(NULL) | |
fb87d707 | 62 | , feopMim(0.9) |
63 | , feopMax(1.3) | |
c2690925 | 64 | { |
65 | // | |
66 | // Constructor | |
67 | // | |
68 | } | |
69 | ||
70 | //___________________________________________________________________ | |
71 | AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c): | |
72 | AliHFEpidBase("") | |
73 | , fPID(NULL) | |
fb87d707 | 74 | , feopMim(0.9) |
75 | , feopMax(1.3) | |
c2690925 | 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 | //___________________________________________________________________ | |
8c1c76e9 | 116 | Bool_t AliHFEpidEMCAL::InitializePID(Int_t /*run*/){ |
c2690925 | 117 | // |
118 | // InitializePID: EMCAL experts have to implement code here | |
119 | // | |
120 | return kTRUE; | |
121 | } | |
122 | ||
123 | //___________________________________________________________________ | |
e156c3bb | 124 | Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const |
c2690925 | 125 | { //Function to a return a code indicating whether or not an electron candidate is selected |
126 | // | |
127 | // | |
fb87d707 | 128 | |
63bffcf1 | 129 | if(track==NULL)return 0; |
130 | ||
8c1c76e9 | 131 | if(!fkPIDResponse) return 0; |
c2690925 | 132 | AliDebug(2, "PID object available"); |
133 | // EMCal not fESDpid (s.s Feb. 11) | |
63bffcf1 | 134 | |
135 | const AliVTrack *trk = dynamic_cast<const AliVTrack *>(track->GetRecTrack()); | |
136 | if (trk == NULL) return 0; | |
e156c3bb | 137 | //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis; |
63bffcf1 | 138 | if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0; |
c2690925 | 139 | AliDebug(2, "Track Has EMCAL PID"); |
e156c3bb | 140 | |
c9f1d8ef | 141 | if(pidqa) |
e156c3bb | 142 | pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID); |
c2690925 | 143 | // not QA for EMCal, will be added (s.s Feb. 11) |
144 | ||
8c1c76e9 | 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 | ||
a8ef1999 | 163 | //printf("eop cuts org; %g; %g \n",feopMim,feopMax); |
164 | //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut); | |
8c1c76e9 | 165 | |
e156c3bb | 166 | Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?) |
c2690925 | 167 | AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop)); |
8c1c76e9 | 168 | AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax)); |
c2690925 | 169 | Int_t pdg = 0; |
8c1c76e9 | 170 | //if(eop>feopMim && eop<feopMax){ |
171 | //if(eop>feopMim && eop<feopMax){ | |
172 | if(eop>feopMimCut && eop<feopMaxCut){ | |
c2690925 | 173 | pdg = 11; |
8c1c76e9 | 174 | //printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop); |
c9f1d8ef | 175 | if(pidqa) |
e156c3bb | 176 | pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID); |
c2690925 | 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 | ||
e156c3bb | 183 | AliDebug(1, Form("eID %g ; %d \n",eop,pdg)); |
8c1c76e9 | 184 | //printf("eID %g ; %d \n",eop,pdg); |
c2690925 | 185 | |
186 | return pdg; | |
187 | } | |
e156c3bb | 188 | |
189 | ||
8c1c76e9 | 190 | //__________________________________________________________________________ |
191 | Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const | |
192 | { | |
57189f04 | 193 | // max nSig cuts for 10d |
194 | ||
8c1c76e9 | 195 | double maxCut = 0.0; |
196 | const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track); | |
197 | if(esdtrack==NULL)return maxCut; | |
57189f04 | 198 | double pt = esdtrack->Pt(); |
8c1c76e9 | 199 | |
200 | if(flageop<0.5) | |
c9f1d8ef | 201 | { //<--- new para for updated non-liniarity |
a8ef1999 | 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}; | |
57189f04 | 204 | double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); |
205 | double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); | |
8c1c76e9 | 206 | maxCut = mean+3.0*sig; |
207 | } | |
208 | else | |
209 | { | |
a8ef1999 | 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}; | |
57189f04 | 212 | double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); |
213 | double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); | |
8c1c76e9 | 214 | maxCut = mean+3.0*sig; |
215 | } | |
a8ef1999 | 216 | |
217 | //printf("eop cuts org; %g; %g ; %g\n",feopMim,feopMax,pt); | |
218 | ||
8c1c76e9 | 219 | return maxCut; |
220 | } | |
221 | ||
222 | Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const | |
223 | { | |
57189f04 | 224 | // mim nsig cuts for 10d |
8c1c76e9 | 225 | double mimCut = 0.0; |
226 | const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track); | |
227 | if(esdtrack==NULL)return mimCut; | |
57189f04 | 228 | double pt = esdtrack->Pt(); |
8c1c76e9 | 229 | |
230 | if(flageop<0.5) // real | |
231 | { | |
232 | //printf("real"); | |
a8ef1999 | 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}; | |
57189f04 | 235 | double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); |
236 | double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); | |
a8ef1999 | 237 | mimCut = mean-3.0*sig; |
8c1c76e9 | 238 | } |
239 | else // MC | |
240 | { | |
241 | //printf("MC"); | |
a8ef1999 | 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}; | |
57189f04 | 244 | double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); |
245 | double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); | |
a8ef1999 | 246 | mimCut = mean-3.0*sig; |
8c1c76e9 | 247 | } |
248 | return mimCut; | |
249 | } | |
250 | ||
e156c3bb | 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); | |
63bffcf1 | 259 | if(esdtrack==NULL)return feop; |
0ecbfc1b | 260 | const AliESDEvent *evt = esdtrack->GetESDEvent(); |
e156c3bb | 261 | |
8c1c76e9 | 262 | Int_t icl = esdtrack->GetEMCALcluster(); |
e156c3bb | 263 | |
264 | AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl); | |
3a7690e5 | 265 | if(cluster==NULL) {return feop;} |
e156c3bb | 266 | if(!cluster->IsEMCAL()) {return feop;} |
267 | ||
268 | matchclsE = cluster->E(); | |
e156c3bb | 269 | |
270 | if(matchclsE<9999.0) feop = matchclsE/esdtrack->P(); | |
271 | ||
272 | return feop; | |
273 | ||
274 | } | |
275 | ||
276 | ||
c2690925 | 277 | /* |
e156c3bb | 278 | //____________________________________________________________________________________ |
279 | Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const | |
c2690925 | 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 |