]>
Commit | Line | Data |
---|---|---|
c43ed0f6 | 1 | //////////////////////////////////////////////////////////////////////////////// |
2 | // // | |
3 | // AliFemtoEventReaderStandard - the reader class for the Alice ESD, AOD // | |
4 | // the model Kinematics information tailored for the Task framework // | |
5 | // Authors: Adam Kisiel Adam.Kisiel@cern.ch // | |
6 | // // | |
7 | //////////////////////////////////////////////////////////////////////////////// | |
8 | #include "AliFemtoEventReaderStandard.h" | |
9 | ||
10 | #include "TFile.h" | |
11 | #include "TTree.h" | |
12 | #include "TList.h" | |
13 | #include "TBits.h" | |
14 | #include "AliESDEvent.h" | |
15 | #include "AliESDtrack.h" | |
16 | #include "AliESDVertex.h" | |
17 | ||
18 | #include "AliFmPhysicalHelixD.h" | |
19 | #include "AliFmThreeVectorF.h" | |
20 | ||
21 | #include "SystemOfUnits.h" | |
22 | ||
23 | #include "AliFemtoEvent.h" | |
24 | ||
25 | #include "TParticle.h" | |
26 | #include "AliFemtoModelHiddenInfo.h" | |
27 | #include "AliFemtoModelGlobalHiddenInfo.h" | |
28 | #include "AliGenHijingEventHeader.h" | |
29 | #include "AliGenCocktailEventHeader.h" | |
30 | ||
31 | #include "AliAODMCHeader.h" | |
32 | #include "AliAODMCParticle.h" | |
33 | ||
34 | #include "AliVertexerTracks.h" | |
35 | ||
36 | ClassImp(AliFemtoEventReaderStandard) | |
37 | ||
38 | #if !(ST_NO_NAMESPACES) | |
39 | using namespace units; | |
40 | #endif | |
41 | ||
42 | using namespace std; | |
43 | //____________________________ | |
44 | AliFemtoEventReaderStandard::AliFemtoEventReaderStandard(): | |
45 | AliFemtoEventReader(), | |
46 | fFileName(" "), | |
47 | fNumberofEvent(0), | |
48 | fCurEvent(0), | |
49 | fCurFile(0), | |
50 | fESDEvent(0x0), | |
51 | fAODEvent(0x0), | |
52 | fStack(0x0), | |
53 | fGenHeader(0x0), | |
54 | fInputType(kUnknown), | |
55 | fUsePhysicsSel(kFALSE), | |
56 | fSelect(0), | |
69c1c8ff | 57 | fTrackCuts(0x0), |
58 | fUseTPCOnly(kFALSE) | |
c43ed0f6 | 59 | { |
60 | //constructor with 0 parameters , look at default settings | |
61 | } | |
62 | ||
63 | //__________________ | |
64 | AliFemtoEventReaderStandard::AliFemtoEventReaderStandard(const AliFemtoEventReaderStandard& aReader): | |
65 | AliFemtoEventReader(aReader), | |
66 | fFileName(" "), | |
67 | fNumberofEvent(0), | |
68 | fCurEvent(0), | |
69 | fCurFile(0), | |
70 | fESDEvent(0x0), | |
71 | fAODEvent(0x0), | |
72 | fStack(0x0), | |
73 | fGenHeader(0x0), | |
74 | fInputType(kUnknown), | |
75 | fUsePhysicsSel(kFALSE), | |
76 | fSelect(0), | |
69c1c8ff | 77 | fTrackCuts(0x0), |
78 | fUseTPCOnly(kFALSE) | |
c43ed0f6 | 79 | { |
80 | // Copy constructor | |
81 | fCurEvent = aReader.fCurEvent; | |
82 | fCurFile = aReader.fCurFile; | |
83 | fESDEvent = new AliESDEvent(); | |
84 | fAODEvent = new AliAODEvent(); | |
85 | fStack = aReader.fStack; | |
86 | fInputType = aReader.fInputType; | |
87 | fUsePhysicsSel = aReader.fUsePhysicsSel; | |
88 | if (fUsePhysicsSel) fSelect = new AliPhysicsSelection(); | |
89 | fTrackCuts = new AliESDtrackCuts(*(aReader.fTrackCuts)); | |
69c1c8ff | 90 | fUseTPCOnly = aReader.fUseTPCOnly; |
c43ed0f6 | 91 | } |
92 | //__________________ | |
93 | AliFemtoEventReaderStandard::~AliFemtoEventReaderStandard() | |
94 | { | |
95 | //Destructor | |
96 | delete fESDEvent; | |
97 | delete fAODEvent; | |
98 | if (fTrackCuts) delete fTrackCuts; | |
99 | } | |
100 | ||
101 | //__________________ | |
102 | AliFemtoEventReaderStandard& AliFemtoEventReaderStandard::operator=(const AliFemtoEventReaderStandard& aReader) | |
103 | { | |
104 | // Assignment operator | |
105 | if (this == &aReader) | |
106 | return *this; | |
107 | ||
108 | fCurEvent = aReader.fCurEvent; | |
109 | fCurFile = aReader.fCurFile; | |
110 | if (fESDEvent) delete fESDEvent; | |
111 | fESDEvent = new AliESDEvent(); | |
112 | if (fAODEvent) delete fAODEvent; | |
113 | fAODEvent = new AliAODEvent(); | |
114 | fStack = aReader.fStack; | |
115 | fGenHeader = aReader.fGenHeader; | |
116 | fInputType = aReader.fInputType; | |
117 | fUsePhysicsSel = aReader.fUsePhysicsSel; | |
118 | if (fUsePhysicsSel) fSelect = new AliPhysicsSelection(); | |
119 | if (fTrackCuts) delete fTrackCuts; | |
120 | fTrackCuts = new AliESDtrackCuts(*(aReader.fTrackCuts)); | |
69c1c8ff | 121 | fUseTPCOnly = aReader.fUseTPCOnly; |
122 | ||
c43ed0f6 | 123 | return *this; |
124 | } | |
125 | //__________________ | |
126 | // Simple report | |
127 | AliFemtoString AliFemtoEventReaderStandard::Report() | |
128 | { | |
129 | AliFemtoString temp = "\n This is the AliFemtoEventReaderStandard\n"; | |
130 | return temp; | |
131 | } | |
132 | //__________________ | |
133 | AliFemtoEvent* AliFemtoEventReaderStandard::ReturnHbtEvent() | |
134 | { | |
135 | // Get the event, read all the relevant information | |
136 | // and fill the AliFemtoEvent class | |
137 | // Returns a valid AliFemtoEvent | |
138 | AliFemtoEvent *hbtEvent = 0; | |
139 | string tFriendFileName; | |
140 | ||
141 | hbtEvent = new AliFemtoEvent; | |
142 | ||
143 | // Get the friend information | |
144 | cout<<"starting to read event "<<fCurEvent<<endl; | |
145 | if ((fInputType == kESD) || (fInputType == kESDKine)) { | |
146 | if(fESDEvent->GetAliESDOld())fESDEvent->CopyFromOldESD(); | |
147 | ||
148 | if (fUsePhysicsSel) { | |
149 | hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fESDEvent)); | |
150 | if (!(fSelect->IsCollisionCandidate(fESDEvent))) | |
151 | printf("Event not a collision candidate\n"); | |
152 | } | |
153 | else | |
154 | hbtEvent->SetIsCollisionCandidate(kTRUE); | |
155 | } | |
156 | else { | |
157 | hbtEvent->SetIsCollisionCandidate(kTRUE); | |
158 | } | |
159 | ||
c43ed0f6 | 160 | double fV1[3]; |
161 | ||
162 | //setting basic things | |
163 | if ((fInputType == kESD) || (fInputType == kESDKine)) { | |
164 | hbtEvent->SetRunNumber(fESDEvent->GetRunNumber()); | |
165 | hbtEvent->SetMagneticField(fESDEvent->GetMagneticField()*kilogauss);//to check if here is ok | |
166 | hbtEvent->SetZDCN1Energy(fESDEvent->GetZDCN1Energy()); | |
167 | hbtEvent->SetZDCP1Energy(fESDEvent->GetZDCP1Energy()); | |
168 | hbtEvent->SetZDCN2Energy(fESDEvent->GetZDCN2Energy()); | |
169 | hbtEvent->SetZDCP2Energy(fESDEvent->GetZDCP2Energy()); | |
170 | hbtEvent->SetZDCEMEnergy(fESDEvent->GetZDCEMEnergy()); | |
171 | hbtEvent->SetZDCParticipants(fESDEvent->GetZDCParticipants()); | |
172 | hbtEvent->SetTriggerMask(fESDEvent->GetTriggerMask()); | |
173 | hbtEvent->SetTriggerCluster(fESDEvent->GetTriggerCluster()); | |
174 | ||
69c1c8ff | 175 | printf("Got event type %i\n", fESDEvent->GetEventType()); |
176 | ||
c43ed0f6 | 177 | //Vertex |
178 | double fVCov[6]; | |
179 | // if (fUseTPCOnly) { | |
180 | // fESDEvent->GetPrimaryVertexTPC()->GetXYZ(fV1); | |
181 | // fESDEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov); | |
182 | // if (!fESDEvent->GetPrimaryVertexTPC()->GetStatus()) | |
183 | // fVCov[4] = -1001.0; | |
184 | // } | |
185 | // else { | |
186 | if (fESDEvent->GetPrimaryVertex()) { | |
187 | fESDEvent->GetPrimaryVertex()->GetXYZ(fV1); | |
188 | fESDEvent->GetPrimaryVertex()->GetCovMatrix(fVCov); | |
189 | ||
190 | if (!fESDEvent->GetPrimaryVertex()->GetStatus()) { | |
191 | // Get the vertex from SPD | |
192 | fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1); | |
193 | fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov); | |
194 | ||
195 | ||
196 | if (!fESDEvent->GetPrimaryVertexSPD()->GetStatus()) | |
197 | fVCov[4] = -1001.0; | |
198 | else { | |
199 | fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1); | |
200 | fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov); | |
201 | } | |
202 | } | |
203 | } | |
204 | else { | |
205 | if (fESDEvent->GetPrimaryVertexSPD()) { | |
206 | fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1); | |
207 | fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov); | |
208 | } | |
209 | } | |
210 | if ((!fESDEvent->GetPrimaryVertex()) && (!fESDEvent->GetPrimaryVertexSPD())) | |
211 | { | |
212 | cout << "No vertex found !!!" << endl; | |
213 | fV1[0] = 10000.0; | |
214 | fV1[1] = 10000.0; | |
215 | fV1[2] = 10000.0; | |
216 | fVCov[4] = -1001.0; | |
217 | } | |
218 | ||
219 | AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]); | |
220 | ||
221 | hbtEvent->SetPrimVertPos(vertex); | |
222 | hbtEvent->SetPrimVertCov(fVCov); | |
223 | } | |
224 | ||
225 | if ((fInputType == kAOD) || (fInputType == kAODKine)) { | |
226 | hbtEvent->SetRunNumber(fAODEvent->GetRunNumber()); | |
227 | hbtEvent->SetMagneticField(fAODEvent->GetMagneticField()*kilogauss);//to check if here is ok | |
228 | hbtEvent->SetZDCN1Energy(fAODEvent->GetZDCN1Energy()); | |
229 | hbtEvent->SetZDCP1Energy(fAODEvent->GetZDCP1Energy()); | |
230 | hbtEvent->SetZDCN2Energy(fAODEvent->GetZDCN2Energy()); | |
231 | hbtEvent->SetZDCP2Energy(fAODEvent->GetZDCP2Energy()); | |
232 | hbtEvent->SetZDCEMEnergy(fAODEvent->GetZDCEMEnergy(0)); | |
233 | hbtEvent->SetZDCParticipants(0); | |
234 | hbtEvent->SetTriggerMask(fAODEvent->GetTriggerMask()); | |
235 | hbtEvent->SetTriggerCluster(fAODEvent->GetTriggerCluster()); | |
236 | ||
237 | // Primary Vertex position | |
238 | fAODEvent->GetPrimaryVertex()->GetPosition(fV1); | |
239 | ||
240 | AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]); | |
241 | hbtEvent->SetPrimVertPos(vertex); | |
242 | } | |
243 | ||
244 | if ((fInputType == kESDKine) || (fInputType == kAODKine)) { | |
245 | Double_t tReactionPlane = 0; | |
246 | ||
247 | AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader); | |
248 | if (!hdh) { | |
249 | AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader); | |
250 | if (cdh) { | |
251 | TList *tGenHeaders = cdh->GetHeaders(); | |
252 | for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) { | |
253 | hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader); | |
254 | if (hdh) break; | |
255 | } | |
256 | } | |
257 | } | |
258 | if (hdh) | |
259 | { | |
260 | tReactionPlane = hdh->ReactionPlaneAngle(); | |
261 | cout << "Got reaction plane " << tReactionPlane << endl; | |
262 | } | |
263 | ||
264 | hbtEvent->SetReactionPlaneAngle(tReactionPlane); | |
265 | } | |
266 | ||
267 | //starting to reading tracks | |
268 | int nofTracks=0; //number of reconstructed tracks in event | |
269 | if ((fInputType == kESD) || (fInputType == kESDKine)) | |
270 | nofTracks = fESDEvent->GetNumberOfTracks(); | |
271 | else if ((fInputType == kAOD) || (fInputType == kAODKine)) | |
272 | nofTracks = fAODEvent->GetNumberOfTracks(); | |
273 | ||
274 | int realnofTracks=0;//number of track which we use ina analysis | |
275 | ||
276 | TClonesArray *mcP = 0; | |
277 | Int_t *motherids=0; | |
278 | ||
279 | if (fInputType == kAODKine) { | |
280 | // Attempt to access MC header | |
281 | AliAODMCHeader *mcH; | |
282 | mcH = (AliAODMCHeader *) fAODEvent->FindListObject(AliAODMCHeader::StdBranchName()); | |
283 | if (!mcH) { | |
284 | cout << "AOD MC information requested, but no header found!" << endl; | |
285 | } | |
286 | ||
287 | mcP = (TClonesArray *) fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()); | |
288 | if (!mcP) { | |
289 | cout << "AOD MC information requested, but no particle array found!" << endl; | |
290 | } | |
291 | ||
292 | hbtEvent->SetReactionPlaneAngle(fAODEvent->GetHeader()->GetQTheta(0)/2.0); | |
293 | ||
294 | if (mcP) { | |
295 | motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()]; | |
296 | for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0; | |
297 | ||
298 | // Read in mother ids | |
299 | AliAODMCParticle *motherpart; | |
300 | for (int ip=0; ip<mcP->GetEntries(); ip++) { | |
301 | motherpart = (AliAODMCParticle *) mcP->At(ip); | |
302 | if (motherpart->GetDaughter(0) > 0) | |
303 | motherids[motherpart->GetDaughter(0)] = ip; | |
304 | if (motherpart->GetDaughter(1) > 0) | |
305 | motherids[motherpart->GetDaughter(1)] = ip; | |
306 | } | |
307 | } | |
308 | } | |
309 | ||
310 | if (fInputType == kESDKine) { | |
311 | motherids = new Int_t[fStack->GetNtrack()]; | |
312 | for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0; | |
313 | ||
314 | // Read in mother ids | |
315 | TParticle *motherpart; | |
316 | for (int ip=0; ip<fStack->GetNtrack(); ip++) { | |
317 | motherpart = fStack->Particle(ip); | |
318 | if (motherpart->GetDaughter(0) > 0) | |
319 | motherids[motherpart->GetDaughter(0)] = ip; | |
320 | if (motherpart->GetDaughter(1) > 0) | |
321 | motherids[motherpart->GetDaughter(1)] = ip; | |
322 | } | |
323 | } | |
324 | ||
325 | for (int i=0;i<nofTracks;i++) | |
326 | { | |
327 | // cout << "Reading track " << i << endl; | |
328 | bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track | |
329 | ||
330 | AliFemtoTrack* trackCopy = new AliFemtoTrack(); | |
331 | ||
332 | if ((fInputType == kESD) || (fInputType == kESDKine)) { | |
69c1c8ff | 333 | |
334 | AliESDtrack *esdtrack = 0x0; | |
335 | if (fUseTPCOnly) { | |
336 | AliESDtrack *mcp = fESDEvent->GetTrack(i); | |
337 | esdtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESDEvent, mcp->GetID()); | |
338 | // printf("Got %p for track %i | ", esdtrack, mcp->GetID()); | |
339 | } | |
340 | else { | |
341 | esdtrack = fESDEvent->GetTrack(i);//getting next track | |
342 | } | |
c43ed0f6 | 343 | |
69c1c8ff | 344 | if (esdtrack && (fTrackCuts->AcceptTrack(esdtrack))) { |
c43ed0f6 | 345 | |
346 | trackCopy->SetCharge((short)esdtrack->GetSign()); | |
347 | ||
348 | //in aliroot we have AliPID | |
349 | //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon | |
350 | //we use only 5 first | |
351 | double esdpid[5]; | |
352 | esdtrack->GetESDpid(esdpid); | |
353 | trackCopy->SetPidProbElectron(esdpid[0]); | |
354 | trackCopy->SetPidProbMuon(esdpid[1]); | |
355 | trackCopy->SetPidProbPion(esdpid[2]); | |
356 | trackCopy->SetPidProbKaon(esdpid[3]); | |
357 | trackCopy->SetPidProbProton(esdpid[4]); | |
358 | ||
359 | double pxyz[3]; | |
360 | double impact[2]; | |
361 | double covimpact[3]; | |
362 | ||
363 | // if (fUseTPCOnly) { | |
364 | // if (!esdtrack->GetTPCInnerParam()) { | |
365 | // cout << "No TPC inner param !" << endl; | |
366 | // delete trackCopy; | |
367 | // continue; | |
368 | // } | |
369 | ||
370 | // AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam()); | |
371 | // param->GetXYZ(rxyz); | |
372 | // param->PropagateToDCA(fESDEvent->GetPrimaryVertexTPC(), (fESDEvent->GetMagneticField()), 10000, impact, covimpact); | |
373 | // param->GetPxPyPz(pxyz);//reading noconstarined momentum | |
374 | ||
375 | // if (fRotateToEventPlane) { | |
376 | // double tPhi = TMath::ATan2(pxyz[1], pxyz[0]); | |
377 | // double tRad = TMath::Hypot(pxyz[0], pxyz[1]); | |
378 | ||
379 | // pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane); | |
380 | // pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane); | |
381 | // } | |
382 | ||
383 | // AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); | |
384 | // if (v.mag() < 0.0001) { | |
385 | // // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl; | |
386 | // delete trackCopy; | |
387 | // continue; | |
388 | // } | |
389 | ||
390 | // trackCopy->SetP(v);//setting momentum | |
391 | // trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); | |
392 | ||
393 | // const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]); | |
394 | // const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]); | |
395 | // //setting helix I do not if it is ok | |
396 | // AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fESDEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); | |
397 | // trackCopy->SetHelix(helix); | |
398 | ||
399 | // //some stuff which could be useful | |
400 | // trackCopy->SetImpactD(impact[0]); | |
401 | // trackCopy->SetImpactZ(impact[1]); | |
402 | // trackCopy->SetCdd(covimpact[0]); | |
403 | // trackCopy->SetCdz(covimpact[1]); | |
404 | // trackCopy->SetCzz(covimpact[2]); | |
405 | // trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact)); | |
406 | ||
407 | // delete param; | |
408 | // } | |
409 | // else { | |
69c1c8ff | 410 | if (fUseTPCOnly) |
411 | tGoodMomentum=esdtrack->GetPxPyPz(pxyz); | |
412 | else | |
413 | tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum | |
414 | // printf("Got good momentum %i\n", tGoodMomentum); | |
c43ed0f6 | 415 | |
416 | AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); | |
69c1c8ff | 417 | if (v.Mag() < 0.0001) { |
c43ed0f6 | 418 | |
419 | delete trackCopy; | |
420 | continue; | |
421 | } | |
422 | ||
423 | trackCopy->SetP(v);//setting momentum | |
424 | trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); | |
425 | const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]); | |
426 | const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]); | |
427 | //setting helix I do not if it is ok | |
428 | AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fESDEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); | |
429 | trackCopy->SetHelix(helix); | |
430 | ||
431 | //some stuff which could be useful | |
432 | float imp[2]; | |
433 | float cim[3]; | |
434 | esdtrack->GetImpactParameters(imp,cim); | |
435 | ||
436 | impact[0] = imp[0]; | |
437 | impact[1] = imp[1]; | |
438 | covimpact[0] = cim[0]; | |
439 | covimpact[1] = cim[1]; | |
440 | covimpact[2] = cim[2]; | |
441 | ||
442 | trackCopy->SetImpactD(impact[0]); | |
443 | trackCopy->SetImpactZ(impact[1]); | |
444 | trackCopy->SetCdd(covimpact[0]); | |
445 | trackCopy->SetCdz(covimpact[1]); | |
446 | trackCopy->SetCzz(covimpact[2]); | |
447 | trackCopy->SetSigmaToVertex(AliESDtrackCuts::GetSigmaToVertex(esdtrack)); | |
448 | ||
449 | trackCopy->SetTrackId(esdtrack->GetID()); | |
450 | trackCopy->SetFlags(esdtrack->GetStatus()); | |
451 | trackCopy->SetLabel(esdtrack->GetLabel()); | |
452 | ||
453 | trackCopy->SetITSchi2(esdtrack->GetITSchi2()); | |
454 | trackCopy->SetITSncls(esdtrack->GetNcls(0)); | |
455 | trackCopy->SetTPCchi2(esdtrack->GetTPCchi2()); | |
456 | trackCopy->SetTPCncls(esdtrack->GetTPCNcls()); | |
457 | trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF()); | |
458 | trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class | |
459 | trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); | |
460 | ||
461 | trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap()); | |
462 | trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap()); | |
463 | ||
464 | double xtpc[3]; | |
465 | esdtrack->GetInnerXYZ(xtpc); | |
466 | xtpc[2] -= fV1[2]; | |
467 | trackCopy->SetNominalTPCEntrancePoint(xtpc); | |
468 | ||
469 | esdtrack->GetOuterXYZ(xtpc); | |
470 | xtpc[2] -= fV1[2]; | |
471 | trackCopy->SetNominalTPCExitPoint(xtpc); | |
472 | ||
473 | int indexes[3]; | |
474 | for (int ik=0; ik<3; ik++) { | |
475 | indexes[ik] = esdtrack->GetKinkIndex(ik); | |
476 | } | |
477 | trackCopy->SetKinkIndexes(indexes); | |
478 | ||
479 | // Fill the hidden information with the simulated data | |
480 | if (fInputType == kESDKine) { | |
481 | if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) { | |
482 | TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel())); | |
483 | ||
484 | // Check the mother information | |
485 | ||
486 | // Using the new way of storing the freeze-out information | |
487 | // Final state particle is stored twice on the stack | |
488 | // one copy (mother) is stored with original freeze-out information | |
489 | // and is not tracked | |
490 | // the other one (daughter) is stored with primary vertex position | |
491 | // and is tracked | |
492 | ||
493 | // Freeze-out coordinates | |
494 | double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; | |
495 | fpx = tPart->Vx() - fV1[0]; | |
496 | fpy = tPart->Vy() - fV1[1]; | |
497 | fpz = tPart->Vz() - fV1[2]; | |
498 | fpt = tPart->T(); | |
499 | ||
500 | AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); | |
501 | tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz); | |
502 | ||
503 | fpx *= 1e13; | |
504 | fpy *= 1e13; | |
505 | fpz *= 1e13; | |
506 | fpt *= 1e13; | |
507 | ||
508 | if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) { | |
509 | TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]); | |
510 | // Check if this is the same particle stored twice on the stack | |
511 | if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) { | |
512 | // It is the same particle | |
513 | // Read in the original freeze-out information | |
514 | // and convert it from to [fm] | |
515 | ||
516 | // EPOS style | |
517 | fpx = mother->Vx()*1e13*0.197327; | |
518 | fpy = mother->Vy()*1e13*0.197327; | |
519 | fpz = mother->Vz()*1e13*0.197327; | |
520 | fpt = mother->T() *1e13*0.197327; | |
521 | ||
522 | ||
523 | // Therminator style | |
524 | // fpx = mother->Vx()*1e13; | |
525 | // fpy = mother->Vy()*1e13; | |
526 | // fpz = mother->Vz()*1e13; | |
527 | // fpt = mother->T() *1e13*3e10; | |
528 | ||
529 | } | |
530 | } | |
531 | ||
532 | tInfo->SetPDGPid(tPart->GetPdgCode()); | |
533 | ||
534 | tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz()); | |
535 | Double_t mass2 = (tPart->Energy() *tPart->Energy() - | |
536 | tPart->Px()*tPart->Px() - | |
537 | tPart->Py()*tPart->Py() - | |
538 | tPart->Pz()*tPart->Pz()); | |
539 | if (mass2>0.0) | |
540 | tInfo->SetMass(TMath::Sqrt(mass2)); | |
541 | else | |
542 | tInfo->SetMass(0.0); | |
543 | ||
544 | tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); | |
545 | trackCopy->SetHiddenInfo(tInfo); | |
546 | } | |
547 | else { | |
548 | AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); | |
549 | tInfo->SetMass(0.0); | |
550 | double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; | |
551 | fpx = fV1[0]*1e13; | |
552 | fpy = fV1[1]*1e13; | |
553 | fpz = fV1[2]*1e13; | |
554 | fpt = 0.0; | |
555 | tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); | |
556 | ||
557 | tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]); | |
558 | ||
559 | trackCopy->SetHiddenInfo(tInfo); | |
560 | } | |
561 | } | |
562 | // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl; | |
563 | } | |
564 | else | |
565 | tGoodMomentum = false; | |
69c1c8ff | 566 | |
567 | if (fUseTPCOnly) | |
568 | if (esdtrack) delete esdtrack; | |
c43ed0f6 | 569 | } |
570 | ||
571 | if ((fInputType == kAOD) || (fInputType == kAODKine)) { | |
572 | // Read in the normal AliAODTracks | |
573 | const AliAODTrack *aodtrack=fAODEvent->GetTrack(i); // getting the AODtrack directly | |
574 | ||
575 | // if (!aodtrack->TestFilterBit(fFilterBit)) | |
576 | // continue; | |
577 | ||
578 | CopyAODtoFemtoTrack(aodtrack, trackCopy); | |
579 | ||
580 | if (mcP) { | |
581 | // Fill the hidden information with the simulated data | |
582 | // Int_t pLabel = aodtrack->GetLabel(); | |
583 | AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel()))); | |
584 | ||
585 | AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); | |
586 | double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; | |
587 | if (!tPart) { | |
588 | fpx = fV1[0]; | |
589 | fpy = fV1[1]; | |
590 | fpz = fV1[2]; | |
591 | tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz); | |
592 | tInfo->SetPDGPid(0); | |
593 | tInfo->SetTrueMomentum(0.0, 0.0, 0.0); | |
594 | tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0); | |
595 | tInfo->SetMass(0); | |
596 | } | |
597 | else { | |
598 | // Check the mother information | |
599 | ||
600 | // Using the new way of storing the freeze-out information | |
601 | // Final state particle is stored twice on the stack | |
602 | // one copy (mother) is stored with original freeze-out information | |
603 | // and is not tracked | |
604 | // the other one (daughter) is stored with primary vertex position | |
605 | // and is tracked | |
606 | ||
607 | // Freeze-out coordinates | |
608 | fpx = tPart->Xv() - fV1[0]; | |
609 | fpy = tPart->Yv() - fV1[1]; | |
610 | fpz = tPart->Zv() - fV1[2]; | |
611 | // fpt = tPart->T(); | |
612 | ||
613 | tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz); | |
614 | ||
615 | fpx *= 1e13; | |
616 | fpy *= 1e13; | |
617 | fpz *= 1e13; | |
618 | // fpt *= 1e13; | |
619 | ||
620 | // cout << "Looking for mother ids " << endl; | |
621 | if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) { | |
622 | // cout << "Got mother id" << endl; | |
623 | AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]); | |
624 | // Check if this is the same particle stored twice on the stack | |
625 | if (mother) { | |
626 | if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) { | |
627 | // It is the same particle | |
628 | // Read in the original freeze-out information | |
629 | // and convert it from to [fm] | |
630 | ||
631 | // EPOS style | |
632 | // fpx = mother->Xv()*1e13*0.197327; | |
633 | // fpy = mother->Yv()*1e13*0.197327; | |
634 | // fpz = mother->Zv()*1e13*0.197327; | |
635 | // fpt = mother->T() *1e13*0.197327*0.5; | |
636 | ||
637 | ||
638 | // Therminator style | |
639 | fpx = mother->Xv()*1e13; | |
640 | fpy = mother->Yv()*1e13; | |
641 | fpz = mother->Zv()*1e13; | |
642 | // fpt = mother->T() *1e13*3e10; | |
643 | ||
644 | } | |
645 | } | |
646 | } | |
647 | ||
648 | // if (fRotateToEventPlane) { | |
649 | // double tPhi = TMath::ATan2(fpy, fpx); | |
650 | // double tRad = TMath::Hypot(fpx, fpy); | |
651 | ||
652 | // fpx = tRad*TMath::Cos(tPhi - tReactionPlane); | |
653 | // fpy = tRad*TMath::Sin(tPhi - tReactionPlane); | |
654 | // } | |
655 | ||
656 | tInfo->SetPDGPid(tPart->GetPdgCode()); | |
657 | ||
658 | // if (fRotateToEventPlane) { | |
659 | // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px()); | |
660 | // double tRad = TMath::Hypot(tPart->Px(), tPart->Py()); | |
661 | ||
662 | // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane), | |
663 | // tRad*TMath::Sin(tPhi - tReactionPlane), | |
664 | // tPart->Pz()); | |
665 | // } | |
666 | // else | |
667 | tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz()); | |
668 | Double_t mass2 = (tPart->E() *tPart->E() - | |
669 | tPart->Px()*tPart->Px() - | |
670 | tPart->Py()*tPart->Py() - | |
671 | tPart->Pz()*tPart->Pz()); | |
672 | if (mass2>0.0) | |
673 | tInfo->SetMass(TMath::Sqrt(mass2)); | |
674 | else | |
675 | tInfo->SetMass(0.0); | |
676 | ||
677 | tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); | |
678 | } | |
679 | trackCopy->SetHiddenInfo(tInfo); | |
680 | } | |
681 | ||
682 | double pxyz[3]; | |
683 | aodtrack->PxPyPz(pxyz);//reading noconstarined momentum | |
684 | const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]); | |
685 | // Check the sanity of the tracks - reject zero momentum tracks | |
69c1c8ff | 686 | if (ktP.Mag() == 0) { |
c43ed0f6 | 687 | delete trackCopy; |
688 | continue; | |
689 | } | |
690 | ||
69c1c8ff | 691 | |
c43ed0f6 | 692 | } |
693 | ||
694 | ||
695 | //decision if we want this track | |
696 | //if we using diffrent labels we want that this label was use for first time | |
697 | //if we use hidden info we want to have match between sim data and ESD | |
698 | if (tGoodMomentum==true) | |
699 | { | |
700 | hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis | |
701 | realnofTracks++;//real number of tracks | |
702 | // delete trackCopy; | |
703 | } | |
704 | else | |
705 | { | |
706 | delete trackCopy; | |
707 | } | |
708 | ||
709 | } | |
710 | ||
711 | if (motherids) | |
b5ed6afb | 712 | delete [] motherids; |
c43ed0f6 | 713 | |
714 | hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event | |
715 | fCurEvent++; | |
716 | // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl; | |
717 | return hbtEvent; | |
718 | } | |
719 | //___________________ | |
720 | void AliFemtoEventReaderStandard::SetESDSource(AliESDEvent *aESD) | |
721 | { | |
722 | // The chain loads the ESD for us | |
723 | // You must provide the address where it can be found | |
724 | fESDEvent = aESD; | |
725 | } | |
726 | //___________________ | |
727 | void AliFemtoEventReaderStandard::SetAODSource(AliAODEvent *aAOD) | |
728 | { | |
729 | // The chain loads the ESD for us | |
730 | // You must provide the address where it can be found | |
731 | fAODEvent = aAOD; | |
732 | } | |
733 | //___________________ | |
734 | void AliFemtoEventReaderStandard::SetStackSource(AliStack *aStack) | |
735 | { | |
736 | // The chain loads the stack for us | |
737 | // You must provide the address where it can be found | |
738 | fStack = aStack; | |
739 | } | |
740 | void AliFemtoEventReaderStandard::SetInputType(AliFemtoInputType aInput) | |
741 | { | |
742 | // Set the proper input type | |
743 | fInputType = aInput; | |
744 | } | |
745 | //___________________ | |
746 | void AliFemtoEventReaderStandard::SetGenEventHeader(AliGenEventHeader *aGenHeader) | |
747 | { | |
748 | // The chain loads the generator event header for us | |
749 | // You must provide the address where it can be found | |
750 | fGenHeader = aGenHeader; | |
751 | } | |
752 | ||
753 | void AliFemtoEventReaderStandard::SetUsePhysicsSelection(const bool usephysics) | |
754 | { | |
755 | fUsePhysicsSel = usephysics; | |
756 | if (!fSelect) fSelect = new AliPhysicsSelection(); | |
757 | } | |
758 | ||
759 | void AliFemtoEventReaderStandard::SetESDTrackCuts(AliESDtrackCuts *esdcuts) | |
760 | { | |
761 | // Set external ESD track cuts | |
762 | fTrackCuts = esdcuts; | |
763 | } | |
764 | ||
69c1c8ff | 765 | void AliFemtoEventReaderStandard::SetUseTPCOnly(const bool usetpconly) |
766 | { | |
767 | // Set flag to use TPC only tracks | |
768 | fUseTPCOnly = usetpconly; | |
769 | } | |
770 | ||
c43ed0f6 | 771 | void AliFemtoEventReaderStandard::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack, |
772 | AliFemtoTrack *tFemtoTrack) | |
773 | { | |
774 | // Copy the track information from the AOD into the internal AliFemtoTrack | |
775 | // If it exists, use the additional information from the PWG2 AOD | |
776 | ||
777 | // Primary Vertex position | |
778 | double fV1[3]; | |
779 | fAODEvent->GetPrimaryVertex()->GetPosition(fV1); | |
780 | ||
781 | tFemtoTrack->SetCharge(tAodTrack->Charge()); | |
782 | ||
783 | //in aliroot we have AliPID | |
784 | //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon | |
785 | //we use only 5 first | |
786 | ||
787 | // AOD pid has 10 components | |
788 | double aodpid[10]; | |
789 | tAodTrack->GetPID(aodpid); | |
790 | tFemtoTrack->SetPidProbElectron(aodpid[0]); | |
791 | tFemtoTrack->SetPidProbMuon(aodpid[1]); | |
792 | tFemtoTrack->SetPidProbPion(aodpid[2]); | |
793 | tFemtoTrack->SetPidProbKaon(aodpid[3]); | |
794 | tFemtoTrack->SetPidProbProton(aodpid[4]); | |
795 | ||
796 | double pxyz[3]; | |
797 | tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum | |
798 | AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); | |
799 | tFemtoTrack->SetP(v);//setting momentum | |
800 | tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); | |
801 | const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]); | |
802 | //setting track helix | |
803 | const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]); | |
804 | AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fAODEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge())); | |
805 | tFemtoTrack->SetHelix(helix); | |
806 | ||
807 | // Flags | |
808 | tFemtoTrack->SetTrackId(tAodTrack->GetID()); | |
809 | tFemtoTrack->SetFlags(1); | |
810 | tFemtoTrack->SetLabel(tAodTrack->GetLabel()); | |
811 | ||
812 | // Track quality information | |
813 | float covmat[6]; | |
814 | tAodTrack->GetCovMatrix(covmat); | |
815 | tFemtoTrack->SetImpactD(covmat[0]); | |
816 | tFemtoTrack->SetImpactZ(covmat[2]); | |
817 | tFemtoTrack->SetCdd(covmat[3]); | |
818 | tFemtoTrack->SetCdz(covmat[4]); | |
819 | tFemtoTrack->SetCzz(covmat[5]); | |
820 | // This information is only available in the ESD | |
821 | // We put in fake values or reasonable estimates | |
822 | tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF()); | |
823 | tFemtoTrack->SetITSncls(1); | |
824 | tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF()); | |
825 | tFemtoTrack->SetTPCncls(1); | |
826 | tFemtoTrack->SetTPCnclsF(1); | |
827 | tFemtoTrack->SetTPCsignalN(1); | |
828 | tFemtoTrack->SetTPCsignalS(1); | |
829 | ||
830 | TBits tAllTrue; | |
831 | TBits tAllFalse; | |
832 | tAllTrue.ResetAllBits(kTRUE); | |
833 | tAllFalse.ResetAllBits(kFALSE); | |
834 | ||
835 | ||
836 | // If not use dummy values | |
837 | tFemtoTrack->SetTPCClusterMap(tAllTrue); | |
838 | tFemtoTrack->SetTPCSharedMap(tAllFalse); | |
839 | ||
840 | double xtpc[3] = {0,0,0}; | |
841 | tFemtoTrack->SetNominalTPCEntrancePoint(xtpc); | |
842 | tFemtoTrack->SetNominalTPCExitPoint(xtpc); | |
843 | ||
844 | int indexes[3]; | |
845 | for (int ik=0; ik<3; ik++) { | |
846 | indexes[ik] = 0; | |
847 | } | |
848 | tFemtoTrack->SetKinkIndexes(indexes); | |
849 | } | |
850 | ||
851 | AliAODMCParticle* AliFemtoEventReaderStandard::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel) | |
852 | { | |
853 | if (aLabel < 0) return 0; | |
854 | AliAODMCParticle *aodP; | |
855 | Int_t posstack = 0; | |
856 | if (aLabel > mcP->GetEntries()) | |
857 | posstack = mcP->GetEntries(); | |
858 | else | |
859 | posstack = aLabel; | |
860 | ||
861 | aodP = (AliAODMCParticle *) mcP->At(posstack); | |
862 | if (aodP->GetLabel() > posstack) { | |
863 | do { | |
864 | aodP = (AliAODMCParticle *) mcP->At(posstack); | |
865 | if (aodP->GetLabel() == aLabel) return aodP; | |
866 | posstack--; | |
867 | } | |
868 | while (posstack > 0); | |
869 | } | |
870 | else { | |
871 | do { | |
872 | aodP = (AliAODMCParticle *) mcP->At(posstack); | |
873 | if (aodP->GetLabel() == aLabel) return aodP; | |
874 | posstack++; | |
875 | } | |
876 | while (posstack < mcP->GetEntries()); | |
877 | } | |
878 | ||
879 | return 0; | |
880 | } |