1 ////////////////////////////////////////////////////////////////////////////////
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 ///
9 ////////////////////////////////////////////////////////////////////////////////
14 *Revision 1.1 2007/05/25 12:42:54 akisiel
15 *Adding a reader for the Kine information
17 *Revision 1.2 2007/05/22 09:01:42 akisiel
18 *Add the possibiloity to save cut settings in the ROOT file
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
23 *Revision 1.5 2007/05/03 09:45:20 akisiel
24 *Fixing Effective C++ warnings
26 *Revision 1.4 2007/04/27 07:28:34 akisiel
27 *Remove event number reading due to interface changes
29 *Revision 1.3 2007/04/27 07:25:16 akisiel
30 *Make revisions needed for compilation from the main AliRoot tree
32 *Revision 1.1.1.1 2007/04/25 15:38:41 panos
33 *Importing the HBT code dir
37 #include "AliFemtoEventReaderESDKine.h"
41 #include "AliESDEvent.h"
42 #include "AliESDtrack.h"
43 #include "AliESDVertex.h"
45 //#include "AliAODParticle.h"
46 #include "TParticle.h"
48 //#include "TSystem.h"
50 #include "AliFmPhysicalHelixD.h"
51 #include "AliFmThreeVectorF.h"
53 #include "SystemOfUnits.h"
55 #include "AliFemtoEvent.h"
58 #include "TParticle.h"
59 #include "AliFemtoModelHiddenInfo.h"
60 #include "AliFemtoModelGlobalHiddenInfo.h"
61 #include "AliGenHijingEventHeader.h"
63 ClassImp(AliFemtoEventReaderESDKine)
65 #if !(ST_NO_NAMESPACES)
66 using namespace units;
70 //____________________________
71 //constructor with 0 parameters , look at default settings
72 AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
85 AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
86 AliFemtoEventReader(),
98 fInputFile = aReader.fInputFile;
99 fFileName = aReader.fFileName;
100 fConstrained = aReader.fConstrained;
101 fNumberofEvent = aReader.fNumberofEvent;
102 fCurEvent = aReader.fCurEvent;
103 fEvent = new AliESDEvent();
107 AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
110 //delete fListOfFiles;
113 if (fRunLoader) delete fRunLoader;
117 AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
119 // assignment operator
120 if (this == &aReader)
123 fInputFile = aReader.fInputFile;
124 fFileName = aReader.fFileName;
125 fConstrained = aReader.fConstrained;
126 fNumberofEvent = aReader.fNumberofEvent;
127 fCurEvent = aReader.fCurEvent;
128 fCurRLEvent = aReader.fCurRLEvent;
129 if (fTree) delete fTree;
130 // fTree = aReader.fTree->CloneTree();
131 if (fEvent) delete fEvent;
132 fEvent = new AliESDEvent();
133 if (fRunLoader) delete fRunLoader;
134 fRunLoader = new AliRunLoader(*aReader.fRunLoader);
139 AliFemtoString AliFemtoEventReaderESDKine::Report()
141 // create reader report
142 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
147 void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
149 //setting the name of file where names of ESD file are written
150 //it takes only this files which have good trees
152 fInputFile=string(inputFile);
153 cout<<"Input File set on "<<fInputFile<<endl;
154 ifstream infile(inputFile);
156 fTree = new TChain("esdTree");
158 if(infile.good()==true)
160 //checking if all give files have good tree inside
161 while (infile.eof()==false)
163 infile.getline(buffer,256);
164 //ifstream test_file(buffer);
165 TFile *esdFile=TFile::Open(buffer,"READ");
168 TTree* tree = (TTree*) esdFile->Get("esdTree");
171 cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
172 fTree->AddFile(buffer);
182 void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
184 fConstrained=constrained;
187 bool AliFemtoEventReaderESDKine::GetConstrained() const
192 AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
194 // read in a next hbt event from the chain
195 // convert it to AliFemtoEvent and return
196 // for further analysis
197 AliFemtoEvent *hbtEvent = 0;
198 TString tGAliceFilename;
200 if (fCurEvent==fNumberofEvent)//open next file
202 if (fNumberofEvent == 0) {
203 fEvent=new AliESDEvent();
206 // fEsdFile=TFile::Open(fFileName.c_str(),"READ");
207 // fTree = (TTree*) fEsdFile->Get("esdTree");
209 fTree->SetBranchStatus("MuonTracks*",0);
210 fTree->SetBranchStatus("PmdTracks*",0);
211 fTree->SetBranchStatus("TrdTracks*",0);
212 fTree->SetBranchStatus("V0s*",0);
213 fTree->SetBranchStatus("Cascades*",0);
214 fTree->SetBranchStatus("Kinks*",0);
215 fTree->SetBranchStatus("CaloClusters*",0);
216 fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
217 fTree->SetBranchStatus("ESDfriend*",0);
218 fEvent->ReadFromTree(fTree);
220 // chain->SetBranchStatus("*",0);
221 // chain->SetBranchStatus("fUniqueID",1);
222 // chain->SetBranchStatus("fTracks",1);
223 // chain->SetBranchStatus("fTracks.*",1);
224 // chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
225 // fTree->SetBranchStatus("fTracks.fCalibContainer",0);
228 fNumberofEvent=fTree->GetEntries();
230 if (fNumberofEvent == 0) {
231 cout<<"no event in input "<<endl;
236 cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
238 // simulation data reading setup
241 else //no more data to read
243 cout<<"no more files "<<hbtEvent<<endl;
248 cout<<"starting to read event "<<fCurEvent<<endl;
249 fTree->GetEvent(fCurEvent);//getting next event
250 // vector<int> tLabelTable;//to check labels
252 cout << "fFileName is " << fFileName.Data() << endl;
253 cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
254 if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
255 fFileName = fTree->GetCurrentFile()->GetName();
256 tGAliceFilename = fFileName;
257 tGAliceFilename.ReplaceAll("AliESDs","galice");
258 cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
259 if (fRunLoader) delete fRunLoader;
260 fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
263 cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
266 if(fRunLoader->LoadHeader())
268 cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
271 fRunLoader->LoadKinematics();
275 fRunLoader->GetEvent(fCurRLEvent);
276 AliStack* tStack = 0x0;
277 tStack = fRunLoader->Stack();
279 hbtEvent = new AliFemtoEvent;
280 //setting basic things
281 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
282 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
283 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
284 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
285 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
286 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
287 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
288 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
289 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
290 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
291 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
292 hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
298 fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
299 fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
300 if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
304 fEvent->GetPrimaryVertex()->GetXYZ(fV1);
305 fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
306 if (!fEvent->GetPrimaryVertex()->GetStatus())
310 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
311 hbtEvent->SetPrimVertPos(vertex);
312 hbtEvent->SetPrimVertCov(fVCov);
314 AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
316 Double_t tReactionPlane = 0;
319 tReactionPlane = hdh->ReactionPlaneAngle();
321 //starting to reading tracks
322 int nofTracks=0; //number of reconstructed tracks in event
323 nofTracks=fEvent->GetNumberOfTracks();
324 int realnofTracks=0;//number of track which we use ina analysis
327 motherids = new Int_t[fStack->GetNtrack()];
328 for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
330 // Read in mother ids
331 TParticle *motherpart;
332 for (int ip=0; ip<fStack->GetNtrack(); ip++) {
333 motherpart = fStack->Particle(ip);
334 if (motherpart->GetDaughter(0) > 0)
335 motherids[motherpart->GetDaughter(0)] = ip;
336 if (motherpart->GetDaughter(1) > 0)
337 motherids[motherpart->GetDaughter(1)] = ip;
339 // if (motherpart->GetPdgCode() == 211) {
340 // cout << "Mother " << ip << " has daughters "
341 // << motherpart->GetDaughter(0) << " "
342 // << motherpart->GetDaughter(1) << " "
343 // << motherpart->Vx() << " "
344 // << motherpart->Vy() << " "
345 // << motherpart->Vz() << " "
351 for (int i=0;i<nofTracks;i++)
353 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
355 AliFemtoTrack* trackCopy = new AliFemtoTrack();
356 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
357 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
359 trackCopy->SetCharge((short)esdtrack->GetSign());
361 //in aliroot we have AliPID
362 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
363 //we use only 5 first
365 esdtrack->GetESDpid(esdpid);
366 trackCopy->SetPidProbElectron(esdpid[0]);
367 trackCopy->SetPidProbMuon(esdpid[1]);
368 trackCopy->SetPidProbPion(esdpid[2]);
369 trackCopy->SetPidProbKaon(esdpid[3]);
370 trackCopy->SetPidProbProton(esdpid[4]);
378 if (!esdtrack->GetTPCInnerParam()) {
383 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
385 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
386 param->GetPxPyPz(pxyz);//reading noconstarined momentum
388 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
389 if (v.mag() < 0.0001) {
390 // cout << "Found 0 momentum ???? " <<endl;
394 trackCopy->SetP(v);//setting momentum
395 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
397 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
398 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
399 //setting helix I do not if it is ok
400 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
401 trackCopy->SetHelix(helix);
403 //some stuff which could be useful
404 trackCopy->SetImpactD(impact[0]);
405 trackCopy->SetImpactZ(impact[1]);
406 trackCopy->SetCdd(covimpact[0]);
407 trackCopy->SetCdz(covimpact[1]);
408 trackCopy->SetCzz(covimpact[2]);
409 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
414 if (fConstrained==true)
415 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
417 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
419 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
420 if (v.mag() < 0.0001) {
421 // cout << "Found 0 momentum ???? " <<endl;
425 trackCopy->SetP(v);//setting momentum
426 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
427 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
428 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
429 //setting helix I do not if it is ok
430 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
431 trackCopy->SetHelix(helix);
433 //some stuff which could be useful
436 esdtrack->GetImpactParameters(imp,cim);
440 covimpact[0] = cim[0];
441 covimpact[1] = cim[1];
442 covimpact[2] = cim[2];
444 trackCopy->SetImpactD(impact[0]);
445 trackCopy->SetImpactZ(impact[1]);
446 trackCopy->SetCdd(covimpact[0]);
447 trackCopy->SetCdz(covimpact[1]);
448 trackCopy->SetCzz(covimpact[2]);
449 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
452 trackCopy->SetTrackId(esdtrack->GetID());
453 trackCopy->SetFlags(esdtrack->GetStatus());
454 trackCopy->SetLabel(esdtrack->GetLabel());
456 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
457 trackCopy->SetITSncls(esdtrack->GetNcls(0));
458 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
459 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
460 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
461 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
462 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
465 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
466 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
469 esdtrack->GetInnerXYZ(xtpc);
471 trackCopy->SetNominalTPCEntrancePoint(xtpc);
473 esdtrack->GetOuterXYZ(xtpc);
475 trackCopy->SetNominalTPCExitPoint(xtpc);
478 for (int ik=0; ik<3; ik++) {
479 indexes[ik] = esdtrack->GetKinkIndex(ik);
481 trackCopy->SetKinkIndexes(indexes);
483 // Fill the hidden information with the simulated data
484 TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
486 // Check the mother information
488 // Using the new way of storing the freeze-out information
489 // Final state particle is stored twice on the stack
490 // one copy (mother) is stored with original freeze-out information
491 // and is not tracked
492 // the other one (daughter) is stored with primary vertex position
495 // Freeze-out coordinates
496 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
497 fpx = tPart->Vx() - fV1[0];
498 fpy = tPart->Vy() - fV1[1];
499 fpz = tPart->Vz() - fV1[2];
502 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
503 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
505 if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
506 TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
507 // Check if this is the same particle stored twice on the stack
508 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
509 // It is the same particle
510 // Read in the original freeze-out information
511 // and convert it from to [fm]
512 fpx = mother->Vx()*1e13;
513 fpy = mother->Vy()*1e13;
514 fpz = mother->Vz()*1e13;
515 fpt = mother->T()*1e13*3e10;
520 tInfo->SetPDGPid(tPart->GetPdgCode());
521 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
522 Double_t mass2 = (tPart->Energy() *tPart->Energy() -
523 tPart->Px()*tPart->Px() -
524 tPart->Py()*tPart->Py() -
525 tPart->Pz()*tPart->Pz());
527 tInfo->SetMass(TMath::Sqrt(mass2));
531 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
532 trackCopy->SetHiddenInfo(tInfo);
534 //decision if we want this track
535 //if we using diffrent labels we want that this label was use for first time
536 //if we use hidden info we want to have match between sim data and ESD
537 if (tGoodMomentum==true)
539 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
540 realnofTracks++;//real number of tracks
550 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
553 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
554 if (fCurEvent== fNumberofEvent)//if end of current file close all
561 //____________________________________________________________________
562 Float_t AliFemtoEventReaderESDKine::GetSigmaToVertex(const AliESDtrack* esdTrack)
564 // Calculates the number of sigma to the vertex.
576 bRes[0] = TMath::Sqrt(bCov[0]);
577 bRes[1] = TMath::Sqrt(bCov[2]);
579 // -----------------------------------
580 // How to get to a n-sigma cut?
582 // The accumulated statistics from 0 to d is
584 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
585 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
587 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
588 // Can this be expressed in a different way?
590 if (bRes[0] == 0 || bRes[1] ==0)
593 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
595 // stupid rounding problem screws up everything:
596 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
597 if (TMath::Exp(-d * d / 2) < 1e-10)
600 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);