Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESDChainKine.cxx
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"
13 #include "TList.h"
14 #include "AliESDEvent.h"
15 #include "AliESDtrack.h"
16 #include "AliESDVertex.h"
17
18 #include "AliMultiplicity.h"
19 #include "AliCentrality.h"
20 #include "AliEventplane.h"
21
22 #include "AliFmPhysicalHelixD.h"
23 #include "AliFmThreeVectorF.h"
24
25 #include "SystemOfUnits.h"
26
27 #include "AliFemtoEvent.h"
28
29 #include "TParticle.h"
30 #include "AliFemtoModelHiddenInfo.h"
31 #include "AliFemtoModelGlobalHiddenInfo.h"
32 #include "AliGenHijingEventHeader.h"
33 #include "AliGenCocktailEventHeader.h"
34
35 #include "AliVertexerTracks.h"
36
37 #include "AliPID.h"
38
39 ClassImp(AliFemtoEventReaderESDChainKine)
40
41 #if !(ST_NO_NAMESPACES)
42   using namespace units;
43 #endif
44
45 using namespace std;
46 //____________________________
47 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
48   fFileName(" "),
49   fConstrained(true),
50   fReadInner(false),
51   fUseTPCOnly(false),
52   fNumberofEvent(0),
53   fCurEvent(0),
54   fCurFile(0),
55   fEvent(0x0),
56   fStack(0x0),
57   fGenHeader(0x0),
58   fTrackType(kGlobal),
59   fEstEventMult(kV0Centrality),
60   fRotateToEventPlane(0),
61   fESDpid(0),
62   fIsPidOwner(0)
63
64 {
65   //constructor with 0 parameters , look at default settings 
66 }
67
68 //__________________
69 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
70   AliFemtoEventReader(aReader),
71   fFileName(" "),
72   fConstrained(true),
73   fReadInner(false),
74   fUseTPCOnly(false),
75   fNumberofEvent(0),
76   fCurEvent(0),
77   fCurFile(0),
78   fEvent(0x0),
79   fStack(0x0),
80   fGenHeader(0x0),
81   fTrackType(kGlobal),
82   fEstEventMult(kV0Centrality),
83   fRotateToEventPlane(0),
84   fESDpid(0),
85   fIsPidOwner(0)
86 {
87   // Copy constructor
88   fConstrained = aReader.fConstrained;
89   fReadInner = aReader.fReadInner;
90   fUseTPCOnly = aReader.fUseTPCOnly;
91   fNumberofEvent = aReader.fNumberofEvent;
92   fCurEvent = aReader.fCurEvent;
93   fCurFile = aReader.fCurFile;
94   fEvent = new AliESDEvent();
95   fStack = aReader.fStack;
96   fTrackType = aReader.fTrackType;
97   fEstEventMult = aReader.fEstEventMult;
98   fRotateToEventPlane = aReader.fRotateToEventPlane;
99   fESDpid = aReader.fESDpid;
100   fIsPidOwner = aReader.fIsPidOwner;
101 }
102 //__________________
103 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
104 {
105   //Destructor
106   delete fEvent;
107 }
108
109 //__________________
110 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
111 {
112   // Assignment operator
113   if (this == &aReader)
114     return *this;
115
116   fConstrained = aReader.fConstrained;
117   fUseTPCOnly = aReader.fUseTPCOnly;
118   fNumberofEvent = aReader.fNumberofEvent;
119   fCurEvent = aReader.fCurEvent;
120   fCurFile = aReader.fCurFile;
121   if (fEvent) delete fEvent;
122   fEvent = new AliESDEvent();
123   fStack = aReader.fStack;
124   fGenHeader = aReader.fGenHeader;
125   fRotateToEventPlane = aReader.fRotateToEventPlane;
126   fESDpid = aReader.fESDpid;
127   fIsPidOwner = aReader.fIsPidOwner;
128
129   return *this;
130 }
131 //__________________
132 // Simple report
133 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
134 {
135   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
136   return temp;
137 }
138
139 //__________________
140 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
141 {
142   // Select whether to read constrained or not constrained momentum
143   fConstrained=constrained;
144 }
145 //__________________
146 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
147 {
148   // Check whether we read constrained or not constrained momentum
149   return fConstrained;
150 }
151 //__________________
152 void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner)
153 {
154   fReadInner=readinner;
155 }
156
157 bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const
158 {
159   return fReadInner;
160 }
161
162 //__________________
163 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
164 {
165   fUseTPCOnly=usetpconly;
166 }
167
168 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
169 {
170   return fUseTPCOnly;
171 }
172 void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType)
173 {
174   fEstEventMult = aType;
175 }
176 //__________________
177 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
178 {
179   // Get the event, read all the relevant information
180   // and fill the AliFemtoEvent class
181   // Returns a valid AliFemtoEvent
182   AliFemtoEvent *hbtEvent = 0;
183   string tFriendFileName;
184
185   // Get the friend information
186   cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
187   //  fEvent->SetESDfriend(fEventFriend);
188         
189   hbtEvent = new AliFemtoEvent;
190   //setting basic things
191   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
192   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
193   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
194   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
195   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
196   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
197   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
198   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
199   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
200   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
201   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
202   //hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
203
204   if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
205       (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
206       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
207       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
208     hbtEvent->SetTriggerCluster(1);
209   else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
210            (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
211     hbtEvent->SetTriggerCluster(2);
212   else 
213     hbtEvent->SetTriggerCluster(0);
214
215   //Vertex
216   double fV1[3];
217   double fVCov[6];
218   if (fUseTPCOnly) {
219     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
220     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
221     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
222       fVCov[4] = -1001.0;
223   }
224   else {
225     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
226     fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
227     if (!fEvent->GetPrimaryVertex()->GetStatus())
228       fVCov[4] = -1001.0;
229   }
230
231   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
232   hbtEvent->SetPrimVertPos(vertex);
233   hbtEvent->SetPrimVertCov(fVCov);
234
235   Double_t tReactionPlane = 0;
236
237   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);  
238   if (!hdh) {
239     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
240     if (cdh) {
241       TList *tGenHeaders = cdh->GetHeaders();
242       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
243         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);     
244         if (hdh) break;
245       }
246     }
247   }
248
249   if (hdh)
250     {
251       tReactionPlane = hdh->ReactionPlaneAngle();
252       cout << "Got reaction plane " << tReactionPlane << endl;
253     }
254
255   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
256
257   Int_t spdetaonecount = 0;
258   
259   for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++) 
260     if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
261       spdetaonecount++;
262
263   //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
264   hbtEvent->SetSPDMult(spdetaonecount);
265
266   //starting to reading tracks
267   int nofTracks=0;  //number of reconstructed tracks in event
268   nofTracks=fEvent->GetNumberOfTracks();
269   int realnofTracks=0;//number of track which we use ina analysis
270
271   int tNormMult = 0;
272   int tNormMultPos = 0;
273   int tNormMultNeg = 0;
274
275   Float_t tTotalPt = 0.0;
276
277   Float_t b[2];
278   Float_t bCov[3];
279
280   Int_t tTracklet=0, tITSTPC=0, tITSPure=0;
281   
282   fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
283   
284   hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
285   hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
286   //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
287   hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
288
289   Int_t *motherids;
290   if (fStack) {
291     motherids = new Int_t[fStack->GetNtrack()];
292     for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
293
294     // Read in mother ids
295     TParticle *motherpart;
296     for (int ip=0; ip<fStack->GetNtrack(); ip++) {
297       motherpart = fStack->Particle(ip);
298       if (motherpart->GetDaughter(0) > 0)
299         motherids[motherpart->GetDaughter(0)] = ip;
300       if (motherpart->GetDaughter(1) > 0)
301         motherids[motherpart->GetDaughter(1)] = ip;
302       
303 //     if (motherpart->GetPdgCode() == 211) {
304 //       cout << "Mother " << ip << " has daughters " 
305 //         << motherpart->GetDaughter(0) << " " 
306 //         << motherpart->GetDaughter(1) << " " 
307 //         << motherpart->Vx() << " " 
308 //         << motherpart->Vy() << " " 
309 //         << motherpart->Vz() << " " 
310 //         << endl;
311       
312 //     }
313     }
314   }
315   else {
316     printf ("No Stack ???\n");
317     delete hbtEvent;
318     return 0;
319   }
320
321   for (int i=0;i<nofTracks;i++)
322     {
323       //      cout << "Reading track " << i << endl;
324       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
325                 
326         
327       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
328       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
329
330
331       if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
332           (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
333         if (esdtrack->GetTPCNcls() > 70) 
334           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
335             if (TMath::Abs(esdtrack->Eta()) < 1.2) {
336               esdtrack->GetImpactParameters(b,bCov);
337               if ((b[0]<0.2) && (b[1] < 0.25)) {
338                 tNormMult++;
339                 tTotalPt += esdtrack->Pt();
340               }
341             }
342           }
343       }
344       else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
345         if (esdtrack->GetTPCNcls() > 100) 
346           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
347             if (TMath::Abs(esdtrack->Eta()) < 1.2) {
348               esdtrack->GetImpactParameters(b,bCov);
349               if ((b[0]<2.4) && (b[1] < 3.2)) {
350                 tNormMult++;
351                 tTotalPt += esdtrack->Pt();
352               }
353             }
354           }
355       }
356       
357       hbtEvent->SetZDCEMEnergy(tTotalPt);
358       //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
359       //        if (esdtrack->GetTPCNcls() > 80) 
360       //          if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0) 
361       //            if (esdtrack->GetConstrainedParam())
362       //              if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
363       //                if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
364       //                  if (esdtrack->GetSign() > 0)
365       //                    tNormMultPos++;
366       //                  else if (esdtrack->GetSign() < 0)
367       //                    tNormMultNeg--;
368       //                }
369
370       // If reading ITS-only tracks, reject all with TPC
371       if (fTrackType == kITSOnly) {
372         if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
373         if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
374         if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
375         UChar_t iclm = esdtrack->GetITSClusterMap();
376         Int_t incls = 0;
377         for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
378         if (incls<=3) {
379           if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
380           continue;
381         }
382       }
383
384
385
386       AliFemtoTrack* trackCopy = new AliFemtoTrack();
387       trackCopy->SetCharge((short)esdtrack->GetSign());
388
389       //in aliroot we have AliPID 
390       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
391       //we use only 5 first
392       double esdpid[5];
393       esdtrack->GetESDpid(esdpid);
394       trackCopy->SetPidProbElectron(esdpid[0]);
395       trackCopy->SetPidProbMuon(esdpid[1]);
396       trackCopy->SetPidProbPion(esdpid[2]);
397       trackCopy->SetPidProbKaon(esdpid[3]);
398       trackCopy->SetPidProbProton(esdpid[4]);
399
400       esdpid[0] = -100000.0;
401       esdpid[1] = -100000.0;
402       esdpid[2] = -100000.0;
403       esdpid[3] = -100000.0;
404       esdpid[4] = -100000.0;
405
406       double tTOF = 0.0;
407
408       if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
409         tTOF = esdtrack->GetTOFsignal();
410         esdtrack->GetIntegratedTimes(esdpid);
411       }
412
413       trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
414
415
416       //////  TPC ////////////////////////////////////////////
417       float nsigmaTPCK=-1000.;                                                        
418       float nsigmaTPCPi=-1000.;                                                        
419       float nsigmaTPCP=-1000.;  
420
421      if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
422         nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
423         nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
424         nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
425
426       }
427       trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
428       trackCopy->SetNSigmaTPCK(nsigmaTPCK);
429       trackCopy->SetNSigmaTPCP(nsigmaTPCP);
430
431       ///// TOF ///////////////////////////////////////////////
432
433       float vp=-1000.;
434       float nsigmaTOFPi=-1000.;
435       float nsigmaTOFK=-1000.;
436       float nsigmaTOFP=-1000.;
437
438         if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
439             (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
440             (esdtrack->GetStatus()&AliESDtrack::kTIME) &&
441             !(esdtrack->GetStatus()&AliESDtrack::kTOFmismatch))
442           {
443
444             //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
445             //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
446             //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
447             // collect info from ESDpid class
448           
449               if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
450               
451               
452               double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
453               
454               nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
455               nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
456               nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
457               
458               Double_t len=esdtrack->GetIntegratedLength();
459               Double_t tof=esdtrack->GetTOFsignal();
460               if(tof > 0.) vp=len/tof/0.03;
461             }
462           }
463
464       trackCopy->SetVTOF(vp);
465       trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
466       trackCopy->SetNSigmaTOFK(nsigmaTOFK);
467       trackCopy->SetNSigmaTOFP(nsigmaTOFP);
468
469                                                 
470       double pxyz[3];
471       double rxyz[3];
472       double impact[2];
473       double covimpact[3];
474       
475       if (fUseTPCOnly) {
476         if (!esdtrack->GetTPCInnerParam()) {
477           cout << "No TPC inner param !" << endl;
478           delete trackCopy;
479           continue;
480         }
481
482         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
483         param->GetXYZ(rxyz);
484         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
485         param->GetPxPyPz(pxyz);//reading noconstarined momentum
486
487           if (fReadInner == true) {
488             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
489             tInfo->SetPDGPid(211);
490             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
491             tInfo->SetMass(0.13957);
492             //    tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
493             //    tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
494             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
495             trackCopy->SetHiddenInfo(tInfo);
496           }
497
498         if (fRotateToEventPlane) {
499           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
500           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
501           
502           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
503           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
504         }
505
506         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
507         if (v.Mag() < 0.0001) {
508           //      cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
509           delete trackCopy;
510           continue;
511         }
512         trackCopy->SetP(v);//setting momentum
513         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
514
515         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
516         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
517         //setting helix I do not if it is ok
518         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
519         trackCopy->SetHelix(helix);
520
521         //some stuff which could be useful 
522         trackCopy->SetImpactD(impact[0]);
523         trackCopy->SetImpactZ(impact[1]);
524         trackCopy->SetCdd(covimpact[0]);
525         trackCopy->SetCdz(covimpact[1]);
526         trackCopy->SetCzz(covimpact[2]);
527         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));       
528
529         delete param;
530       }
531       else {
532           if (fReadInner == true) {
533           
534             if (esdtrack->GetTPCInnerParam()) {
535               AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
536               param->GetXYZ(rxyz);
537               //            param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
538               param->GetPxPyPz(pxyz);//reading noconstarined momentum
539               delete param;
540             
541               AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
542               tInfo->SetPDGPid(211);
543               tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
544               tInfo->SetMass(0.13957);
545               //            tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
546               //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
547               tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
548               trackCopy->SetHiddenInfo(tInfo);
549             }
550           }
551
552
553         if(fTrackType == kGlobal) {
554         if (fConstrained==true)             
555           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
556         else
557           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
558         }
559           else if (fTrackType == kTPCOnly) {
560             if (esdtrack->GetTPCInnerParam())
561               esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
562             else {
563               delete trackCopy;
564               continue;
565             }
566           }
567           else if (fTrackType == kITSOnly) {
568             if (fConstrained==true)                 
569               tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
570             else
571               tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
572           }
573
574         if (fRotateToEventPlane) {
575           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
576           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
577           
578           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
579           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
580         }
581
582         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
583         if (v.Mag() < 0.0001) {
584           //      cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
585           delete trackCopy;
586           continue;
587         }
588
589         trackCopy->SetP(v);//setting momentum
590         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
591         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
592         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
593         //setting helix I do not if it is ok
594         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
595         trackCopy->SetHelix(helix);
596         
597         //some stuff which could be useful 
598         float imp[2];
599         float cim[3];
600         esdtrack->GetImpactParameters(imp,cim);
601
602         impact[0] = imp[0];
603         impact[1] = imp[1];
604         covimpact[0] = cim[0];
605         covimpact[1] = cim[1];
606         covimpact[2] = cim[2];
607
608         trackCopy->SetImpactD(impact[0]);
609         trackCopy->SetImpactZ(impact[1]);
610         trackCopy->SetCdd(covimpact[0]);
611         trackCopy->SetCdz(covimpact[1]);
612         trackCopy->SetCzz(covimpact[2]);
613         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
614       }
615                 
616       trackCopy->SetTrackId(esdtrack->GetID());
617       trackCopy->SetFlags(esdtrack->GetStatus());
618       trackCopy->SetLabel(esdtrack->GetLabel());
619                 
620       trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
621       trackCopy->SetITSncls(esdtrack->GetNcls(0));     
622       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
623       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
624       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
625       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
626       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
627
628
629       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
630       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
631
632       double xtpc[3];
633       esdtrack->GetInnerXYZ(xtpc);
634       xtpc[2] -= fV1[2];
635       trackCopy->SetNominalTPCEntrancePoint(xtpc);
636
637       esdtrack->GetOuterXYZ(xtpc);
638       xtpc[2] -= fV1[2];
639       trackCopy->SetNominalTPCExitPoint(xtpc);
640
641       int indexes[3];
642       for (int ik=0; ik<3; ik++) {
643         indexes[ik] = esdtrack->GetKinkIndex(ik);
644       }
645       trackCopy->SetKinkIndexes(indexes);
646
647
648       // Fill the hidden information with the simulated data
649       if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
650         TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
651
652         // Check the mother information
653         
654         // Using the new way of storing the freeze-out information
655         // Final state particle is stored twice on the stack
656         // one copy (mother) is stored with original freeze-out information
657         //   and is not tracked
658         // the other one (daughter) is stored with primary vertex position
659         //   and is tracked
660         
661         // Freeze-out coordinates
662         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
663         fpx = tPart->Vx() - fV1[0];
664         fpy = tPart->Vy() - fV1[1];
665         fpz = tPart->Vz() - fV1[2];
666         fpt = tPart->T();
667         
668         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
669         tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
670         
671         fpx *= 1e13;
672         fpy *= 1e13;
673         fpz *= 1e13;
674         fpt *= 1e13;
675         
676         //      cout << "Looking for mother ids " << endl;
677         if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
678           //    cout << "Got mother id" << endl;
679           TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
680           // Check if this is the same particle stored twice on the stack
681           if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
682             // It is the same particle
683             // Read in the original freeze-out information
684             // and convert it from to [fm]
685             
686             // EPOS style 
687             fpx = mother->Vx()*1e13*0.197327;
688             fpy = mother->Vy()*1e13*0.197327;
689             fpz = mother->Vz()*1e13*0.197327;
690             fpt = mother->T() *1e13*0.197327;
691             
692             
693             // Therminator style 
694 //          fpx = mother->Vx()*1e13;
695 //          fpy = mother->Vy()*1e13;
696 //          fpz = mother->Vz()*1e13;
697 //          fpt = mother->T() *1e13*3e10;
698             
699           }
700         }
701         
702         if (fRotateToEventPlane) {
703           double tPhi = TMath::ATan2(fpy, fpx);
704           double tRad = TMath::Hypot(fpx, fpy);
705           
706           fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
707           fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
708         }
709         
710         tInfo->SetPDGPid(tPart->GetPdgCode());
711         
712         if (fRotateToEventPlane) {
713           double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
714           double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
715           
716           tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
717                                  tRad*TMath::Sin(tPhi - tReactionPlane),
718                                  tPart->Pz());
719         }
720         else
721           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
722         Double_t mass2 = (tPart->Energy() *tPart->Energy() -
723                           tPart->Px()*tPart->Px() -
724                           tPart->Py()*tPart->Py() -
725                           tPart->Pz()*tPart->Pz());
726         if (mass2>0.0)
727           tInfo->SetMass(TMath::Sqrt(mass2));
728         else 
729           tInfo->SetMass(0.0);
730         
731         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
732         trackCopy->SetHiddenInfo(tInfo);
733       }
734       else {
735         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
736         tInfo->SetMass(0.0);
737         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
738         fpx = fV1[0]*1e13;
739         fpy = fV1[1]*1e13;
740         fpz = fV1[2]*1e13;
741         fpt = 0.0;
742         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
743
744         tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
745         
746         trackCopy->SetHiddenInfo(tInfo);
747       }
748       //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl; 
749       
750       //decision if we want this track
751       //if we using diffrent labels we want that this label was use for first time 
752       //if we use hidden info we want to have match between sim data and ESD
753       if (tGoodMomentum==true)
754         {
755           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
756           realnofTracks++;//real number of tracks
757           //      delete trackCopy;
758         }
759       else
760         {
761           delete  trackCopy;
762         }
763                 
764     }
765
766   delete [] motherids;
767   
768   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event   
769
770   AliCentrality *cent = fEvent->GetCentrality();
771   if (cent) {
772     hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
773     //    hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
774     hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
775     //    hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
776
777     if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
778   }
779
780   if (fEstEventMult == kGlobalCount) 
781     hbtEvent->SetNormalizedMult(tNormMult);
782   else if (fEstEventMult == kTracklet)
783     hbtEvent->SetNormalizedMult(tTracklet);
784   else if (fEstEventMult == kITSTPC)
785     hbtEvent->SetNormalizedMult(tITSTPC);
786   else if (fEstEventMult == kITSPure)
787     hbtEvent->SetNormalizedMult(tITSPure);
788   else if (fEstEventMult == kSPDLayer1)
789     hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
790   else if (fEstEventMult == kV0Centrality) {
791     // centrality between 0 (central) and 1 (very peripheral)
792
793     if (cent) {
794       if (cent->GetCentralityPercentile("V0M") < 0.00001)
795         hbtEvent->SetNormalizedMult(-1);
796       else
797         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
798       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
799                              10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
800     }
801   }
802
803   if (tNormMultPos > tNormMultNeg)
804     hbtEvent->SetZDCParticipants(tNormMultPos);
805   else
806     hbtEvent->SetZDCParticipants(tNormMultNeg);
807   
808   AliEventplane* ep = fEvent->GetEventplane();
809   if (ep) {
810     hbtEvent->SetEP(ep);
811     hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
812   }
813
814
815   fCurEvent++;  
816   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
817   return hbtEvent; 
818 }
819 //___________________
820 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
821 {
822   // The chain loads the ESD for us
823   // You must provide the address where it can be found
824   fEvent = aESD;
825 }
826 //___________________
827 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
828 {
829   // The chain loads the stack for us
830   // You must provide the address where it can be found
831   fStack = aStack;
832 }
833 //___________________
834 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
835 {
836   // The chain loads the generator event header for us
837   // You must provide the address where it can be found
838   fGenHeader = aGenHeader;
839 }
840 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
841 {
842   fRotateToEventPlane=dorotate;
843 }
844 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
845 {
846   // Calculates the number of sigma to the vertex.
847
848   Float_t b[2];
849   Float_t bRes[2];
850   Float_t bCov[3];
851
852   b[0] = impact[0];
853   b[1] = impact[1];
854   bCov[0] = covar[0];
855   bCov[1] = covar[1];
856   bCov[2] = covar[2];
857
858   bRes[0] = TMath::Sqrt(bCov[0]);
859   bRes[1] = TMath::Sqrt(bCov[2]);
860
861   // -----------------------------------
862   // How to get to a n-sigma cut?
863   //
864   // The accumulated statistics from 0 to d is
865   //
866   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
867   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
868   //
869   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
870   // Can this be expressed in a different way?
871
872   if (bRes[0] == 0 || bRes[1] ==0)
873     return -1;
874
875   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
876
877   // stupid rounding problem screws up everything:
878   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
879   if (TMath::Exp(-d * d / 2) < 1e-10)
880     return 1000;
881
882   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
883   return d;
884 }
885
886 void AliFemtoEventReaderESDChainKine::SetReadTrackType(ReadTrackType aType)
887 {
888   fTrackType = aType;
889 }
890