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 }