Bring AliFemto up to date with latest code developements
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoEventReaderESDKine.cxx
1 ////////////////////////////////////////////////////////////////////////////////
2 ///                                                                          ///
3 /// AliFemtoEventReaderESDKine - the reader class for the Alice ESD              ///
4 /// Reads in ESD information and converts it into internal AliFemtoEvent     ///
5 /// Reads in AliESDfriend to create shared hit/quality information           ///
6 /// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl                        ///
7 ///          Adam Kisiel kisiel@mps.ohio-state.edu                           ///
8 ///                                                                          ///
9 ////////////////////////////////////////////////////////////////////////////////
10
11 /*
12  *$Id$
13  *$Log$
14  *Revision 1.1  2007/05/25 12:42:54  akisiel
15  *Adding a reader for the Kine information
16  *
17  *Revision 1.2  2007/05/22 09:01:42  akisiel
18  *Add the possibiloity to save cut settings in the ROOT file
19  *
20  *Revision 1.1  2007/05/16 10:22:11  akisiel
21  *Making the directory structure of AliFemto flat. All files go into one common directory
22  *
23  *Revision 1.5  2007/05/03 09:45:20  akisiel
24  *Fixing Effective C++ warnings
25  *
26  *Revision 1.4  2007/04/27 07:28:34  akisiel
27  *Remove event number reading due to interface changes
28  *
29  *Revision 1.3  2007/04/27 07:25:16  akisiel
30  *Make revisions needed for compilation from the main AliRoot tree
31  *
32  *Revision 1.1.1.1  2007/04/25 15:38:41  panos
33  *Importing the HBT code dir
34  *
35  */
36
37 #include "AliFemtoEventReaderESDKine.h"
38
39 #include "TFile.h"
40 #include "TChain.h"
41 #include "AliESDEvent.h"
42 #include "AliESDtrack.h"
43 #include "AliESDVertex.h"
44 #include "AliStack.h"
45 //#include "AliAODParticle.h"
46 #include "TParticle.h"
47
48 //#include "TSystem.h"
49
50 #include "AliFmPhysicalHelixD.h"
51 #include "AliFmThreeVectorF.h"
52
53 #include "SystemOfUnits.h"
54
55 #include "AliFemtoEvent.h"
56
57 #include "TMath.h"
58
59 ClassImp(AliFemtoEventReaderESDKine)
60
61 #if !(ST_NO_NAMESPACES)
62   using namespace units;
63 #endif
64
65 using namespace std;
66 //____________________________
67 //constructor with 0 parameters , look at default settings 
68 AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
69   fInputFile(" "),
70   fFileName(" "),
71   fConstrained(true),
72   fNumberofEvent(0),
73   fCurEvent(0),
74   fCurRLEvent(0),
75   fTree(0x0),
76   fEvent(0x0),
77   fRunLoader(0x0)
78 {
79 }
80
81 AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
82   fInputFile(" "),
83   fFileName(" "),
84   fConstrained(true),
85   fNumberofEvent(0),
86   fCurEvent(0),
87   fCurRLEvent(0),
88   fTree(0x0),
89   fEvent(0x0),
90   fRunLoader(0x0)
91 {
92   // copy constructor
93   fInputFile = aReader.fInputFile;
94   fFileName  = aReader.fFileName;
95   fConstrained = aReader.fConstrained;
96   fNumberofEvent = aReader.fNumberofEvent;
97   fCurEvent = aReader.fCurEvent;
98   fEvent = new AliESDEvent();
99 }
100 //__________________
101 //Destructor
102 AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
103 {
104   // destructor
105   //delete fListOfFiles;
106   delete fTree;
107   delete fEvent;
108   if (fRunLoader) delete fRunLoader;
109 }
110
111 //__________________
112 AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
113 {
114   // assignment operator
115   if (this == &aReader)
116     return *this;
117
118   fInputFile = aReader.fInputFile;
119   fFileName  = aReader.fFileName;
120   fConstrained = aReader.fConstrained;
121   fNumberofEvent = aReader.fNumberofEvent;
122   fCurEvent = aReader.fCurEvent;
123   fCurRLEvent = aReader.fCurRLEvent;
124   if (fTree) delete fTree;
125   //  fTree = aReader.fTree->CloneTree();
126   if (fEvent) delete fEvent;
127   fEvent = new AliESDEvent();
128   if (fRunLoader) delete fRunLoader;
129   fRunLoader = new AliRunLoader(*aReader.fRunLoader);
130
131   return *this;
132 }
133 //__________________
134 AliFemtoString AliFemtoEventReaderESDKine::Report() const
135 {
136   // create reader report
137   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
138   return temp;
139 }
140
141 //__________________
142 void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
143 {
144   //setting the name of file where names of ESD file are written 
145   //it takes only this files which have good trees
146   char buffer[256];
147   fInputFile=string(inputFile);
148   cout<<"Input File set on "<<fInputFile<<endl;
149   ifstream infile(inputFile);
150
151   fTree = new TChain("esdTree");
152
153   if(infile.good()==true)
154     { 
155       //checking if all give files have good tree inside
156       while (infile.eof()==false)
157         {
158           infile.getline(buffer,256);
159           //ifstream test_file(buffer);
160           TFile *esdFile=TFile::Open(buffer,"READ");
161           if (esdFile!=0x0)
162             {   
163               TTree* tree = (TTree*) esdFile->Get("esdTree");
164               if (tree!=0x0)
165                 {
166                   cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
167                   fTree->AddFile(buffer);
168                   delete tree;
169                 }
170               esdFile->Close(); 
171             }
172           delete esdFile;
173         }
174     }
175 }
176
177 void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
178 {
179   fConstrained=constrained;
180 }
181
182 bool AliFemtoEventReaderESDKine::GetConstrained() const
183 {
184   return fConstrained;
185 }
186
187 AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
188 {
189   // read in a next hbt event from the chain
190   // convert it to AliFemtoEvent and return
191   // for further analysis
192   AliFemtoEvent *hbtEvent = 0;
193   TString tGAliceFilename;
194
195   if (fCurEvent==fNumberofEvent)//open next file  
196     {
197       if (fNumberofEvent == 0) {
198         fEvent=new AliESDEvent();
199                 
200           //ESD data
201 //        fEsdFile=TFile::Open(fFileName.c_str(),"READ");
202 //        fTree = (TTree*) fEsdFile->Get("esdTree");                    
203
204           fTree->SetBranchStatus("MuonTracks*",0);
205           fTree->SetBranchStatus("PmdTracks*",0);
206           fTree->SetBranchStatus("TrdTracks*",0);
207           fTree->SetBranchStatus("V0s*",0);
208           fTree->SetBranchStatus("Cascades*",0);
209           fTree->SetBranchStatus("Kinks*",0);
210           fTree->SetBranchStatus("CaloClusters*",0);
211           fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
212           fTree->SetBranchStatus("ESDfriend*",0);
213           fEvent->ReadFromTree(fTree);
214
215 //        chain->SetBranchStatus("*",0);
216 //        chain->SetBranchStatus("fUniqueID",1);
217 //        chain->SetBranchStatus("fTracks",1);
218 //        chain->SetBranchStatus("fTracks.*",1);
219 //        chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
220 //        fTree->SetBranchStatus("fTracks.fCalibContainer",0);
221
222
223         fNumberofEvent=fTree->GetEntries();
224
225         if (fNumberofEvent == 0) {
226           cout<<"no event in input "<<endl;
227           fReaderStatus=1;
228           return hbtEvent; 
229         }
230
231         cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
232         fCurEvent=0;
233         // simulation data reading setup
234         
235       }
236       else //no more data to read
237         {
238           cout<<"no more files "<<hbtEvent<<endl;
239           fReaderStatus=1;
240           return hbtEvent; 
241         }
242     }           
243   cout<<"starting to read event "<<fCurEvent<<endl;
244   fTree->GetEvent(fCurEvent);//getting next event
245   //  vector<int> tLabelTable;//to check labels
246   
247   cout << "fFileName is " << fFileName.Data() << endl;
248   cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
249   if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
250     fFileName = fTree->GetCurrentFile()->GetName();
251     tGAliceFilename = fFileName;
252     tGAliceFilename.ReplaceAll("AliESDs","galice");
253     cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
254     if (fRunLoader) delete fRunLoader;
255     fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
256     if (fRunLoader==0x0)
257       {
258         cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
259         exit(0);
260       }
261     if(fRunLoader->LoadHeader())
262       {
263         cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
264         exit(0);
265       }
266     fRunLoader->LoadKinematics();
267     fCurRLEvent = 0;
268   }
269
270   fRunLoader->GetEvent(fCurRLEvent);
271   AliStack* tStack = 0x0;
272   tStack = fRunLoader->Stack();
273         
274   hbtEvent = new AliFemtoEvent;
275   //setting basic things
276   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
277   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
278   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
279   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
280   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
281   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
282   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
283   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
284   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
285   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
286   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
287   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
288         
289   //Vertex
290   double fV1[3];
291   fEvent->GetVertex()->GetXYZ(fV1);
292
293   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
294   hbtEvent->SetPrimVertPos(vertex);
295         
296   //starting to reading tracks
297   int nofTracks=0;  //number of reconstructed tracks in event
298   nofTracks=fEvent->GetNumberOfTracks();
299   int realnofTracks=0;//number of track which we use ina analysis
300
301   for (int i=0;i<nofTracks;i++)
302     {
303       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
304                 
305       AliFemtoTrack* trackCopy = new AliFemtoTrack();   
306       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
307       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
308
309       trackCopy->SetCharge((short)esdtrack->GetSign());
310
311       //in aliroot we have AliPID 
312       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
313       //we use only 5 first
314       double esdpid[5];
315       esdtrack->GetESDpid(esdpid);
316       trackCopy->SetPidProbElectron(esdpid[0]);
317       trackCopy->SetPidProbMuon(esdpid[1]);
318       trackCopy->SetPidProbPion(esdpid[2]);
319       trackCopy->SetPidProbKaon(esdpid[3]);
320       trackCopy->SetPidProbProton(esdpid[4]);
321                                                 
322       double pxyz[3];
323       if (fConstrained==true)               
324         tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
325       else
326         tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
327       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
328       trackCopy->SetP(v);//setting momentum
329       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
330       const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
331       if (ktP.mag() == 0) {
332         delete trackCopy;
333         continue;
334       }
335       const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
336       //setting helix I do not if it is ok
337       AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
338       trackCopy->SetHelix(helix);
339                 
340       trackCopy->SetTrackId(esdtrack->GetID());
341       trackCopy->SetFlags(esdtrack->GetStatus());
342       trackCopy->SetLabel(esdtrack->GetLabel());
343                 
344       //some stuff which could be useful 
345       float impact[2];
346       float covimpact[3];
347       esdtrack->GetImpactParameters(impact,covimpact);
348       trackCopy->SetImpactD(impact[0]);
349       trackCopy->SetImpactZ(impact[1]);
350       trackCopy->SetCdd(covimpact[0]);
351       trackCopy->SetCdz(covimpact[1]);
352       trackCopy->SetCzz(covimpact[2]);
353       trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
354       trackCopy->SetITSncls(esdtrack->GetNcls(0));     
355       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
356       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
357       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
358       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
359       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
360       trackCopy->SetSigmaToVertex(GetSigmaToVertex(esdtrack));
361
362       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
363       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
364
365       double pvrt[3];
366       fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
367
368       double xtpc[3];
369       esdtrack->GetInnerXYZ(xtpc);
370       xtpc[2] -= pvrt[2];
371       trackCopy->SetNominalTPCEntrancePoint(xtpc);
372
373       esdtrack->GetOuterXYZ(xtpc);
374       xtpc[2] -= pvrt[2];
375       trackCopy->SetNominalTPCExitPoint(xtpc);
376
377       int indexes[3];
378       for (int ik=0; ik<3; ik++) {
379         indexes[ik] = esdtrack->GetKinkIndex(ik);
380       }
381       trackCopy->SetKinkIndexes(indexes);
382
383       // Fill the hidden information with the simulated data
384       TParticle *tPart = tStack->Particle(TMath::Abs(esdtrack->GetLabel()));
385       //      AliAODParticle* tParticle= new AliAODParticle(*tPart,i);
386       AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
387       tInfo->SetPDGPid(tPart->GetPdgCode());
388       tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
389       Double_t mass2 = (tPart->Energy() *tPart->Energy() -
390                         tPart->Px()*tPart->Px() -
391                         tPart->Py()*tPart->Py() -
392                         tPart->Pz()*tPart->Pz());
393       if (mass2>0.0)
394         tInfo->SetMass(TMath::Sqrt(mass2));
395       else 
396         tInfo->SetMass(0.0);
397       tInfo->SetEmissionPoint(tPart->Vx()*1e13, tPart->Vy()*1e13, tPart->Vz()*1e13, tPart->T()*1e13*300000);
398       trackCopy->SetHiddenInfo(tInfo);
399
400       //decision if we want this track
401       //if we using diffrent labels we want that this label was use for first time 
402       //if we use hidden info we want to have match between sim data and ESD
403       if (tGoodMomentum==true)
404         {
405           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
406           realnofTracks++;//real number of tracks
407         }
408       else
409         {
410           delete  trackCopy;
411         }
412                 
413     }
414
415   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event   
416   fCurEvent++;  
417   fCurRLEvent++;
418   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
419   if (fCurEvent== fNumberofEvent)//if end of current file close all
420     {   
421       fTree->Reset(); 
422       delete fTree;
423     }
424   return hbtEvent; 
425 }
426 //____________________________________________________________________
427 Float_t AliFemtoEventReaderESDKine::GetSigmaToVertex(const AliESDtrack* esdTrack)
428 {
429   // Calculates the number of sigma to the vertex.
430
431   Float_t b[2];
432   Float_t bRes[2];
433   Float_t bCov[3];
434   esdTrack->GetImpactParameters(b,bCov);
435 //   if (bCov[0]<=0 || bCov[2]<=0) {
436 //     AliDebug(1, "Estimated b resolution lower or equal zero!");
437 //     bCov[0]=0; bCov[2]=0;
438 //   }
439   bRes[0] = TMath::Sqrt(bCov[0]);
440   bRes[1] = TMath::Sqrt(bCov[2]);
441
442   // -----------------------------------
443   // How to get to a n-sigma cut?
444   //
445   // The accumulated statistics from 0 to d is
446   //
447   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
448   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
449   //
450   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
451   // Can this be expressed in a different way?
452
453   if (bRes[0] == 0 || bRes[1] ==0)
454     return -1;
455
456   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
457
458   // stupid rounding problem screws up everything:
459   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
460   if (TMath::Exp(-d * d / 2) < 1e-10)
461     return 1000;
462
463   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
464   return d;
465 }