Port of changes from v4-07-Release and additional rule conformance
[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 ClassImp(AliFemtoEventReaderESDKine)
58
59 #if !(ST_NO_NAMESPACES)
60   using namespace units;
61 #endif
62
63 using namespace std;
64 //____________________________
65 //constructor with 0 parameters , look at default settings 
66 AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
67   fInputFile(" "),
68   fFileName(" "),
69   fConstrained(true),
70   fNumberofEvent(0),
71   fCurEvent(0),
72   fCurRLEvent(0),
73   fTree(0x0),
74   fEvent(0x0),
75   fRunLoader(0x0)
76 {
77 }
78
79 AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
80   fInputFile(" "),
81   fFileName(" "),
82   fConstrained(true),
83   fNumberofEvent(0),
84   fCurEvent(0),
85   fCurRLEvent(0),
86   fTree(0x0),
87   fEvent(0x0),
88   fRunLoader(0x0)
89 {
90   // copy constructor
91   fInputFile = aReader.fInputFile;
92   fFileName  = aReader.fFileName;
93   fConstrained = aReader.fConstrained;
94   fNumberofEvent = aReader.fNumberofEvent;
95   fCurEvent = aReader.fCurEvent;
96   fEvent = new AliESDEvent();
97 }
98 //__________________
99 //Destructor
100 AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
101 {
102   // destructor
103   //delete fListOfFiles;
104   delete fTree;
105   delete fEvent;
106   if (fRunLoader) delete fRunLoader;
107 }
108
109 //__________________
110 AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
111 {
112   // assignment operator
113   if (this == &aReader)
114     return *this;
115
116   fInputFile = aReader.fInputFile;
117   fFileName  = aReader.fFileName;
118   fConstrained = aReader.fConstrained;
119   fNumberofEvent = aReader.fNumberofEvent;
120   fCurEvent = aReader.fCurEvent;
121   fCurRLEvent = aReader.fCurRLEvent;
122   if (fTree) delete fTree;
123   //  fTree = aReader.fTree->CloneTree();
124   if (fEvent) delete fEvent;
125   fEvent = new AliESDEvent();
126   if (fRunLoader) delete fRunLoader;
127   fRunLoader = new AliRunLoader(*aReader.fRunLoader);
128
129   return *this;
130 }
131 //__________________
132 AliFemtoString AliFemtoEventReaderESDKine::Report() const
133 {
134   // create reader report
135   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
136   return temp;
137 }
138
139 //__________________
140 void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
141 {
142   //setting the name of file where names of ESD file are written 
143   //it takes only this files which have good trees
144   char buffer[256];
145   fInputFile=string(inputFile);
146   cout<<"Input File set on "<<fInputFile<<endl;
147   ifstream infile(inputFile);
148
149   fTree = new TChain("esdTree");
150
151   if(infile.good()==true)
152     { 
153       //checking if all give files have good tree inside
154       while (infile.eof()==false)
155         {
156           infile.getline(buffer,256);
157           //ifstream test_file(buffer);
158           TFile *esdFile=TFile::Open(buffer,"READ");
159           if (esdFile!=0x0)
160             {   
161               TTree* tree = (TTree*) esdFile->Get("esdTree");
162               if (tree!=0x0)
163                 {
164                   cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
165                   fTree->AddFile(buffer);
166                   delete tree;
167                 }
168               esdFile->Close(); 
169             }
170           delete esdFile;
171         }
172     }
173 }
174
175 void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
176 {
177   fConstrained=constrained;
178 }
179
180 bool AliFemtoEventReaderESDKine::GetConstrained() const
181 {
182   return fConstrained;
183 }
184
185 AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
186 {
187   // read in a next hbt event from the chain
188   // convert it to AliFemtoEvent and return
189   // for further analysis
190   AliFemtoEvent *hbtEvent = 0;
191   TString tGAliceFilename;
192
193   if (fCurEvent==fNumberofEvent)//open next file  
194     {
195       if (fNumberofEvent == 0) {
196         fEvent=new AliESDEvent();
197                 
198           //ESD data
199 //        fEsdFile=TFile::Open(fFileName.c_str(),"READ");
200 //        fTree = (TTree*) fEsdFile->Get("esdTree");                    
201
202           fTree->SetBranchStatus("MuonTracks*",0);
203           fTree->SetBranchStatus("PmdTracks*",0);
204           fTree->SetBranchStatus("TrdTracks*",0);
205           fTree->SetBranchStatus("V0s*",0);
206           fTree->SetBranchStatus("Cascades*",0);
207           fTree->SetBranchStatus("Kinks*",0);
208           fTree->SetBranchStatus("CaloClusters*",0);
209           fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
210           fTree->SetBranchStatus("ESDfriend*",0);
211           fEvent->ReadFromTree(fTree);
212
213 //        chain->SetBranchStatus("*",0);
214 //        chain->SetBranchStatus("fUniqueID",1);
215 //        chain->SetBranchStatus("fTracks",1);
216 //        chain->SetBranchStatus("fTracks.*",1);
217 //        chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
218 //        fTree->SetBranchStatus("fTracks.fCalibContainer",0);
219
220
221         fNumberofEvent=fTree->GetEntries();
222
223         if (fNumberofEvent == 0) {
224           cout<<"no event in input "<<endl;
225           fReaderStatus=1;
226           return hbtEvent; 
227         }
228
229         cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
230         fCurEvent=0;
231         // simulation data reading setup
232         
233       }
234       else //no more data to read
235         {
236           cout<<"no more files "<<hbtEvent<<endl;
237           fReaderStatus=1;
238           return hbtEvent; 
239         }
240     }           
241   cout<<"starting to read event "<<fCurEvent<<endl;
242   fTree->GetEvent(fCurEvent);//getting next event
243   //  vector<int> tLabelTable;//to check labels
244   
245   cout << "fFileName is " << fFileName.Data() << endl;
246   cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
247   if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
248     fFileName = fTree->GetCurrentFile()->GetName();
249     tGAliceFilename = fFileName;
250     tGAliceFilename.ReplaceAll("AliESDs","galice");
251     cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
252     if (fRunLoader) delete fRunLoader;
253     fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
254     if (fRunLoader==0x0)
255       {
256         cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
257         exit(0);
258       }
259     if(fRunLoader->LoadHeader())
260       {
261         cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
262         exit(0);
263       }
264     fRunLoader->LoadKinematics();
265     fCurRLEvent = 0;
266   }
267
268   fRunLoader->GetEvent(fCurRLEvent);
269   AliStack* tStack = 0x0;
270   tStack = fRunLoader->Stack();
271         
272   hbtEvent = new AliFemtoEvent;
273   //setting basic things
274   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
275   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
276   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
277   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
278   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
279   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
280   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
281   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
282   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
283   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
284   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
285   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
286         
287   //Vertex
288   double fV1[3];
289   fEvent->GetVertex()->GetXYZ(fV1);
290
291   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
292   hbtEvent->SetPrimVertPos(vertex);
293         
294   //starting to reading tracks
295   int nofTracks=0;  //number of reconstructed tracks in event
296   nofTracks=fEvent->GetNumberOfTracks();
297   int realnofTracks=0;//number of track which we use ina analysis
298
299   for (int i=0;i<nofTracks;i++)
300     {
301       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
302                 
303       AliFemtoTrack* trackCopy = new AliFemtoTrack();   
304       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
305       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
306
307       trackCopy->SetCharge((short)esdtrack->GetSign());
308
309       //in aliroot we have AliPID 
310       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
311       //we use only 5 first
312       double esdpid[5];
313       esdtrack->GetESDpid(esdpid);
314       trackCopy->SetPidProbElectron(esdpid[0]);
315       trackCopy->SetPidProbMuon(esdpid[1]);
316       trackCopy->SetPidProbPion(esdpid[2]);
317       trackCopy->SetPidProbKaon(esdpid[3]);
318       trackCopy->SetPidProbProton(esdpid[4]);
319                                                 
320       double pxyz[3];
321       if (fConstrained==true)               
322         tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
323       else
324         tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
325       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
326       trackCopy->SetP(v);//setting momentum
327       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
328       const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
329       if (ktP.mag() == 0) {
330         delete trackCopy;
331         continue;
332       }
333       const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
334       //setting helix I do not if it is ok
335       AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
336       trackCopy->SetHelix(helix);
337                 
338       trackCopy->SetTrackId(esdtrack->GetID());
339       trackCopy->SetFlags(esdtrack->GetStatus());
340       trackCopy->SetLabel(esdtrack->GetLabel());
341                 
342       //some stuff which could be useful 
343       float impact[2];
344       float covimpact[3];
345       esdtrack->GetImpactParameters(impact,covimpact);
346       trackCopy->SetImpactD(impact[0]);
347       trackCopy->SetImpactZ(impact[1]);
348       trackCopy->SetCdd(covimpact[0]);
349       trackCopy->SetCdz(covimpact[1]);
350       trackCopy->SetCzz(covimpact[2]);
351       trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
352       trackCopy->SetITSncls(esdtrack->GetNcls(0));     
353       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
354       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
355       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
356       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
357       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
358
359       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
360       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
361
362       // Fill the hidden information with the simulated data
363       TParticle *tPart = tStack->Particle(TMath::Abs(esdtrack->GetLabel()));
364       AliAODParticle* tParticle= new AliAODParticle(*tPart,i);
365       AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
366       tInfo->SetPDGPid(tParticle->GetMostProbable());
367       tInfo->SetTrueMomentum(tParticle->Px(), tParticle->Py(), tParticle->Pz());
368       Double_t mass2 = (tParticle->E()*tParticle->E() -
369                          tParticle->Px()*tParticle->Px() -
370                          tParticle->Py()*tParticle->Py() -
371                          tParticle->Pz()*tParticle->Pz());
372       if (mass2>0.0)
373         tInfo->SetMass(TMath::Sqrt(mass2));
374       else 
375         tInfo->SetMass(0.0);
376       trackCopy->SetHiddenInfo(tInfo);
377
378       //decision if we want this track
379       //if we using diffrent labels we want that this label was use for first time 
380       //if we use hidden info we want to have match between sim data and ESD
381       if (tGoodMomentum==true)
382         {
383           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
384           realnofTracks++;//real number of tracks
385         }
386       else
387         {
388           delete  trackCopy;
389         }
390                 
391     }
392
393   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event   
394   fCurEvent++;  
395   fCurRLEvent++;
396   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
397   if (fCurEvent== fNumberofEvent)//if end of current file close all
398     {   
399       fTree->Reset(); 
400       delete fTree;
401     }
402   return hbtEvent; 
403 }