Add possibility to rotate event
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESDChainKine.cxx
CommitLineData
0b3bd1ac 1////////////////////////////////////////////////////////////////////////////////
2// //
3// AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and //
4// the model Kinematics information tailored for the Task framework and the //
5// Reads in AliESDfriend to create shared hit/quality information //
6// Authors: Adam Kisiel kisiel@mps.ohio-state.edu //
7// //
8////////////////////////////////////////////////////////////////////////////////
9#include "AliFemtoEventReaderESDChainKine.h"
10
11#include "TFile.h"
12#include "TTree.h"
f3523544 13#include "TList.h"
0b3bd1ac 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
0b3bd1ac 25#include "TParticle.h"
26#include "AliFemtoModelHiddenInfo.h"
4c399116 27#include "AliFemtoModelGlobalHiddenInfo.h"
aad44575 28#include "AliGenHijingEventHeader.h"
f3523544 29#include "AliGenCocktailEventHeader.h"
0b3bd1ac 30
d6ce11d4 31#include "AliVertexerTracks.h"
32
0b3bd1ac 33ClassImp(AliFemtoEventReaderESDChainKine)
34
35#if !(ST_NO_NAMESPACES)
36 using namespace units;
37#endif
38
39using namespace std;
40//____________________________
41AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
42 fFileName(" "),
43 fConstrained(true),
613e9805 44 fUseTPCOnly(false),
0b3bd1ac 45 fNumberofEvent(0),
46 fCurEvent(0),
47 fCurFile(0),
48 fEvent(0x0),
aad44575 49 fStack(0x0),
f3523544 50 fGenHeader(0x0),
51 fRotateToEventPlane(0)
0b3bd1ac 52{
53 //constructor with 0 parameters , look at default settings
54}
55
56//__________________
57AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
fcda1d4e 58 AliFemtoEventReader(aReader),
0b3bd1ac 59 fFileName(" "),
60 fConstrained(true),
613e9805 61 fUseTPCOnly(false),
0b3bd1ac 62 fNumberofEvent(0),
63 fCurEvent(0),
64 fCurFile(0),
65 fEvent(0x0),
aad44575 66 fStack(0x0),
f3523544 67 fGenHeader(0x0),
68 fRotateToEventPlane(0)
0b3bd1ac 69{
70 // Copy constructor
71 fConstrained = aReader.fConstrained;
613e9805 72 fUseTPCOnly = aReader.fUseTPCOnly;
0b3bd1ac 73 fNumberofEvent = aReader.fNumberofEvent;
74 fCurEvent = aReader.fCurEvent;
75 fCurFile = aReader.fCurFile;
76 fEvent = new AliESDEvent();
77 fStack = aReader.fStack;
f3523544 78 fRotateToEventPlane = aReader.fRotateToEventPlane;
0b3bd1ac 79}
80//__________________
81AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
82{
83 //Destructor
84 delete fEvent;
85}
86
87//__________________
88AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
89{
90 // Assignment operator
91 if (this == &aReader)
92 return *this;
93
94 fConstrained = aReader.fConstrained;
613e9805 95 fUseTPCOnly = aReader.fUseTPCOnly;
0b3bd1ac 96 fNumberofEvent = aReader.fNumberofEvent;
97 fCurEvent = aReader.fCurEvent;
98 fCurFile = aReader.fCurFile;
99 if (fEvent) delete fEvent;
100 fEvent = new AliESDEvent();
101 fStack = aReader.fStack;
aad44575 102 fGenHeader = aReader.fGenHeader;
f3523544 103 fRotateToEventPlane = aReader.fRotateToEventPlane;
0b3bd1ac 104 return *this;
105}
106//__________________
107// Simple report
108AliFemtoString AliFemtoEventReaderESDChainKine::Report()
109{
110 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
111 return temp;
112}
113
114//__________________
115void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
116{
117 // Select whether to read constrained or not constrained momentum
118 fConstrained=constrained;
119}
120//__________________
121bool AliFemtoEventReaderESDChainKine::GetConstrained() const
122{
123 // Check whether we read constrained or not constrained momentum
124 return fConstrained;
125}
126//__________________
127AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
128{
129 // Get the event, read all the relevant information
130 // and fill the AliFemtoEvent class
131 // Returns a valid AliFemtoEvent
132 AliFemtoEvent *hbtEvent = 0;
133 string tFriendFileName;
134
135 // Get the friend information
e7155bf6 136 cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
0b3bd1ac 137 // fEvent->SetESDfriend(fEventFriend);
138
139 hbtEvent = new AliFemtoEvent;
140 //setting basic things
141 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
142 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
143 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
144 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
145 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
146 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
147 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
148 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
149 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
150 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
151 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
152 hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
153
154 //Vertex
155 double fV1[3];
63a5982a 156 double fVCov[6];
613e9805 157 if (fUseTPCOnly) {
158 fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
63a5982a 159 fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
160 if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
161 fVCov[4] = -1001.0;
613e9805 162 }
163 else {
6fc46bc2 164 if (fEvent->GetPrimaryVertex()) {
165 fEvent->GetPrimaryVertex()->GetXYZ(fV1);
166 fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
167
168 if (!fEvent->GetPrimaryVertex()->GetStatus()) {
169 // Get the vertex from SPD
170 fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
171 fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
172
173
174 if (!fEvent->GetPrimaryVertexSPD()->GetStatus())
175 fVCov[4] = -1001.0;
176 else {
177 fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
178 fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
179 }
180 }
181 }
182 else {
183 if (fEvent->GetPrimaryVertexSPD()) {
d6ce11d4 184 fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
185 fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
186 }
187 }
6fc46bc2 188 if ((!fEvent->GetPrimaryVertex()) && (!fEvent->GetPrimaryVertexSPD()))
189 {
190 cout << "No vertex found !!!" << endl;
191 fV1[0] = 10000.0;
192 fV1[1] = 10000.0;
193 fV1[2] = 10000.0;
194 fVCov[4] = -1001.0;
195 }
613e9805 196 }
0b3bd1ac 197
198 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
63a5982a 199
0b3bd1ac 200 hbtEvent->SetPrimVertPos(vertex);
63a5982a 201 hbtEvent->SetPrimVertCov(fVCov);
aad44575 202
aad44575 203 Double_t tReactionPlane = 0;
f3523544 204
205 AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
206 if (!hdh) {
207 AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
208 if (cdh) {
209 TList *tGenHeaders = cdh->GetHeaders();
210 for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
211 hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
212 if (hdh) break;
213 }
214 }
215 }
216
aad44575 217 if (hdh)
218 {
219 tReactionPlane = hdh->ReactionPlaneAngle();
f3523544 220 cout << "Got reaction plane " << tReactionPlane << endl;
aad44575 221 }
0b3bd1ac 222 //starting to reading tracks
223 int nofTracks=0; //number of reconstructed tracks in event
224 nofTracks=fEvent->GetNumberOfTracks();
225 int realnofTracks=0;//number of track which we use ina analysis
226
613e9805 227 Int_t *motherids;
228 motherids = new Int_t[fStack->GetNtrack()];
229 for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
230
231 // Read in mother ids
232 TParticle *motherpart;
233 for (int ip=0; ip<fStack->GetNtrack(); ip++) {
234 motherpart = fStack->Particle(ip);
235 if (motherpart->GetDaughter(0) > 0)
236 motherids[motherpart->GetDaughter(0)] = ip;
237 if (motherpart->GetDaughter(1) > 0)
238 motherids[motherpart->GetDaughter(1)] = ip;
239
240// if (motherpart->GetPdgCode() == 211) {
241// cout << "Mother " << ip << " has daughters "
242// << motherpart->GetDaughter(0) << " "
243// << motherpart->GetDaughter(1) << " "
244// << motherpart->Vx() << " "
245// << motherpart->Vy() << " "
246// << motherpart->Vz() << " "
247// << endl;
248
249// }
250 }
251
0b3bd1ac 252 for (int i=0;i<nofTracks;i++)
253 {
f3523544 254 // cout << "Reading track " << i << endl;
0b3bd1ac 255 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
256
257 AliFemtoTrack* trackCopy = new AliFemtoTrack();
258 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
259 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
260
261 trackCopy->SetCharge((short)esdtrack->GetSign());
262
263 //in aliroot we have AliPID
264 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
265 //we use only 5 first
266 double esdpid[5];
267 esdtrack->GetESDpid(esdpid);
268 trackCopy->SetPidProbElectron(esdpid[0]);
269 trackCopy->SetPidProbMuon(esdpid[1]);
270 trackCopy->SetPidProbPion(esdpid[2]);
271 trackCopy->SetPidProbKaon(esdpid[3]);
272 trackCopy->SetPidProbProton(esdpid[4]);
273
274 double pxyz[3];
613e9805 275 double rxyz[3];
276 double impact[2];
277 double covimpact[3];
278
279 if (fUseTPCOnly) {
280 if (!esdtrack->GetTPCInnerParam()) {
e7155bf6 281 cout << "No TPC inner param !" << endl;
613e9805 282 delete trackCopy;
283 continue;
284 }
285
286 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
287 param->GetXYZ(rxyz);
288 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
289 param->GetPxPyPz(pxyz);//reading noconstarined momentum
0b3bd1ac 290
f3523544 291 if (fRotateToEventPlane) {
292 double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
293 double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
294
295 pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
296 pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
297 }
298
613e9805 299 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
300 if (v.mag() < 0.0001) {
e7155bf6 301 // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
613e9805 302 delete trackCopy;
303 continue;
304 }
f3523544 305
613e9805 306 trackCopy->SetP(v);//setting momentum
307 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
308
309 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
310 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
311 //setting helix I do not if it is ok
312 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
313 trackCopy->SetHelix(helix);
314
315 //some stuff which could be useful
316 trackCopy->SetImpactD(impact[0]);
317 trackCopy->SetImpactZ(impact[1]);
318 trackCopy->SetCdd(covimpact[0]);
319 trackCopy->SetCdz(covimpact[1]);
320 trackCopy->SetCzz(covimpact[2]);
321 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
322
323 delete param;
324 }
325 else {
326 if (fConstrained==true)
327 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
328 else
329 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
330
f3523544 331
332 if (fRotateToEventPlane) {
333 double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
334 double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
335
336 pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
337 pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
338 }
339
613e9805 340 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
341 if (v.mag() < 0.0001) {
e7155bf6 342 // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
613e9805 343 delete trackCopy;
344 continue;
345 }
f3523544 346
613e9805 347 trackCopy->SetP(v);//setting momentum
348 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
349 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
350 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
351 //setting helix I do not if it is ok
352 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
353 trackCopy->SetHelix(helix);
354
355 //some stuff which could be useful
356 float imp[2];
357 float cim[3];
358 esdtrack->GetImpactParameters(imp,cim);
359
360 impact[0] = imp[0];
361 impact[1] = imp[1];
362 covimpact[0] = cim[0];
363 covimpact[1] = cim[1];
364 covimpact[2] = cim[2];
365
366 trackCopy->SetImpactD(impact[0]);
367 trackCopy->SetImpactZ(impact[1]);
368 trackCopy->SetCdd(covimpact[0]);
369 trackCopy->SetCdz(covimpact[1]);
370 trackCopy->SetCzz(covimpact[2]);
371 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
0b3bd1ac 372 }
0b3bd1ac 373
374 trackCopy->SetTrackId(esdtrack->GetID());
375 trackCopy->SetFlags(esdtrack->GetStatus());
376 trackCopy->SetLabel(esdtrack->GetLabel());
377
0b3bd1ac 378 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
379 trackCopy->SetITSncls(esdtrack->GetNcls(0));
380 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
381 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
382 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
383 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
384 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
385
386
387 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
388 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
389
0b3bd1ac 390 double xtpc[3];
391 esdtrack->GetInnerXYZ(xtpc);
613e9805 392 xtpc[2] -= fV1[2];
0b3bd1ac 393 trackCopy->SetNominalTPCEntrancePoint(xtpc);
394
395 esdtrack->GetOuterXYZ(xtpc);
613e9805 396 xtpc[2] -= fV1[2];
0b3bd1ac 397 trackCopy->SetNominalTPCExitPoint(xtpc);
398
399 int indexes[3];
400 for (int ik=0; ik<3; ik++) {
401 indexes[ik] = esdtrack->GetKinkIndex(ik);
402 }
403 trackCopy->SetKinkIndexes(indexes);
404
405 // Fill the hidden information with the simulated data
406 TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
613e9805 407
408 // Check the mother information
409
410 // Using the new way of storing the freeze-out information
411 // Final state particle is stored twice on the stack
412 // one copy (mother) is stored with original freeze-out information
413 // and is not tracked
414 // the other one (daughter) is stored with primary vertex position
415 // and is tracked
416
417 // Freeze-out coordinates
418 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
1f445b96 419 fpx = tPart->Vx() - fV1[0];
420 fpy = tPart->Vy() - fV1[1];
421 fpz = tPart->Vz() - fV1[2];
613e9805 422 fpt = tPart->T();
423
4c399116 424 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
425 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
426
d6ce11d4 427 fpx *= 1e13;
428 fpy *= 1e13;
429 fpz *= 1e13;
f3523544 430 fpt *= 1e13;
d6ce11d4 431
f3523544 432 // cout << "Looking for mother ids " << endl;
613e9805 433 if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
f3523544 434 // cout << "Got mother id" << endl;
613e9805 435 TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
436 // Check if this is the same particle stored twice on the stack
437 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
438 // It is the same particle
439 // Read in the original freeze-out information
440 // and convert it from to [fm]
f3523544 441
442 // EPOS style
443// fpx = mother->Vx()*1e13*0.197327;
444// fpy = mother->Vy()*1e13*0.197327;
445// fpz = mother->Vz()*1e13*0.197327;
446// fpt = mother->T() *1e13*0.197327*0.5;
447
448
449 // Therminator style
613e9805 450 fpx = mother->Vx()*1e13;
451 fpy = mother->Vy()*1e13;
452 fpz = mother->Vz()*1e13;
f3523544 453 fpt = mother->T() *1e13*3e10;
613e9805 454
455 }
456 }
457
f3523544 458 if (fRotateToEventPlane) {
459 double tPhi = TMath::ATan2(fpy, fpx);
460 double tRad = TMath::Hypot(fpx, fpy);
461
462 fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
463 fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
464 }
465
aad44575 466 tInfo->SetPDGPid(tPart->GetPdgCode());
f3523544 467
468 if (fRotateToEventPlane) {
469 double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
470 double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
471
472 tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
473 tRad*TMath::Sin(tPhi - tReactionPlane),
474 tPart->Pz());
475 }
476 else
477 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
aad44575 478 Double_t mass2 = (tPart->Energy() *tPart->Energy() -
479 tPart->Px()*tPart->Px() -
480 tPart->Py()*tPart->Py() -
481 tPart->Pz()*tPart->Pz());
0b3bd1ac 482 if (mass2>0.0)
483 tInfo->SetMass(TMath::Sqrt(mass2));
484 else
485 tInfo->SetMass(0.0);
aad44575 486
613e9805 487 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
0b3bd1ac 488 trackCopy->SetHiddenInfo(tInfo);
f3523544 489
490 // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl;
613e9805 491
0b3bd1ac 492 //decision if we want this track
493 //if we using diffrent labels we want that this label was use for first time
494 //if we use hidden info we want to have match between sim data and ESD
495 if (tGoodMomentum==true)
496 {
497 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
498 realnofTracks++;//real number of tracks
499 // delete trackCopy;
500 }
501 else
502 {
503 delete trackCopy;
504 }
505
506 }
f3523544 507
508 delete motherids;
613e9805 509
0b3bd1ac 510 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
511 fCurEvent++;
e7155bf6 512 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
0b3bd1ac 513 return hbtEvent;
514}
515//___________________
516void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
517{
518 // The chain loads the ESD for us
519 // You must provide the address where it can be found
520 fEvent = aESD;
521}
522//___________________
523void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
524{
525 // The chain loads the stack for us
526 // You must provide the address where it can be found
527 fStack = aStack;
528}
aad44575 529//___________________
530void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
531{
532 // The chain loads the generator event header for us
533 // You must provide the address where it can be found
534 fGenHeader = aGenHeader;
535}
0b3bd1ac 536
613e9805 537//__________________
538void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
539{
540 fUseTPCOnly=usetpconly;
541}
f3523544 542void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
543{
544 fRotateToEventPlane=dorotate;
545}
0b3bd1ac 546
613e9805 547bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
548{
549 return fUseTPCOnly;
550}
551Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
552{
553 // Calculates the number of sigma to the vertex.
554
555 Float_t b[2];
556 Float_t bRes[2];
557 Float_t bCov[3];
558
559 b[0] = impact[0];
560 b[1] = impact[1];
561 bCov[0] = covar[0];
562 bCov[1] = covar[1];
563 bCov[2] = covar[2];
564
565 bRes[0] = TMath::Sqrt(bCov[0]);
566 bRes[1] = TMath::Sqrt(bCov[2]);
567
568 // -----------------------------------
569 // How to get to a n-sigma cut?
570 //
571 // The accumulated statistics from 0 to d is
572 //
573 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
574 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
575 //
576 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
577 // Can this be expressed in a different way?
578
579 if (bRes[0] == 0 || bRes[1] ==0)
580 return -1;
581
582 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
583
584 // stupid rounding problem screws up everything:
585 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
586 if (TMath::Exp(-d * d / 2) < 1e-10)
587 return 1000;
588
589 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
590 return d;
591}