PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderStandard.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoEventReaderStandard - the reader class for the Alice ESD, AOD      //
   4 // the model Kinematics information tailored for the Task framework           //
   5 // Authors: Adam Kisiel Adam.Kisiel@cern.ch                                   //
   6 //                                                                            //
   7 ////////////////////////////////////////////////////////////////////////////////
   8 #include "AliFemtoEventReaderStandard.h"
   9
  10 #include "TFile.h"
  11 #include "TTree.h"
  12 #include "TList.h"
  13 #include "TBits.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDtrack.h"
  16 #include "AliESDVertex.h"
  17
  18 #include "AliFmPhysicalHelixD.h"
  19 #include "AliFmThreeVectorF.h"
  20
  21 #include "SystemOfUnits.h"
  22
  23 #include "AliFemtoEvent.h"
  24
  25 #include "TParticle.h"
  26 #include "AliFemtoModelHiddenInfo.h"
  27 #include "AliFemtoModelGlobalHiddenInfo.h"
  28 #include "AliGenHijingEventHeader.h"
  29 #include "AliGenCocktailEventHeader.h"
  30
  31 #include "AliAODMCHeader.h"
  32 #include "AliAODMCParticle.h"
  33
  34 #include "AliVertexerTracks.h"
  35
  36 ClassImp(AliFemtoEventReaderStandard)
  37
  38 #if !(ST_NO_NAMESPACES)
  39   using namespace units;
  40 #endif
  41
  42 using namespace std;
  43 //____________________________
  44 AliFemtoEventReaderStandard::AliFemtoEventReaderStandard():
  45   AliFemtoEventReader(),
  46   fFileName(" "),
  47   fNumberofEvent(0),
  48   fCurEvent(0),
  49   fCurFile(0),
  50   fESDEvent(0x0),
  51   fAODEvent(0x0),
  52   fStack(0x0),
  53   fGenHeader(0x0),
  54   fInputType(kUnknown),
  55   fUsePhysicsSel(kFALSE),
  56   fSelect(0),
  57   fTrackCuts(0x0),
  58   fUseTPCOnly(kFALSE)
  59 {
  60   //constructor with 0 parameters , look at default settings
  61 }
  62
  63 //__________________
  64 AliFemtoEventReaderStandard::AliFemtoEventReaderStandard(const AliFemtoEventReaderStandard& aReader):
  65   AliFemtoEventReader(aReader),
  66   fFileName(" "),
  67   fNumberofEvent(0),
  68   fCurEvent(0),
  69   fCurFile(0),
  70   fESDEvent(0x0),
  71   fAODEvent(0x0),
  72   fStack(0x0),
  73   fGenHeader(0x0),
  74   fInputType(kUnknown),
  75   fUsePhysicsSel(kFALSE),
  76   fSelect(0),
  77   fTrackCuts(0x0),
  78   fUseTPCOnly(kFALSE)
  79 {
  80   // Copy constructor
  81   fCurEvent = aReader.fCurEvent;
  82   fCurFile = aReader.fCurFile;
  83   fESDEvent = new AliESDEvent();
  84   fAODEvent = new AliAODEvent();
  85   fStack = aReader.fStack;
  86   fInputType = aReader.fInputType;
  87   fUsePhysicsSel = aReader.fUsePhysicsSel;
  88   if (fUsePhysicsSel) fSelect = new AliPhysicsSelection();
  89   fTrackCuts = new AliESDtrackCuts(*(aReader.fTrackCuts));
  90   fUseTPCOnly = aReader.fUseTPCOnly;
  91 }
  92 //__________________
  93 AliFemtoEventReaderStandard::~AliFemtoEventReaderStandard()
  94 {
  95   //Destructor
  96   delete fESDEvent;
  97   delete fAODEvent;
  98   if (fTrackCuts) delete fTrackCuts;
  99 }
 100
 101 //__________________
 102 AliFemtoEventReaderStandard& AliFemtoEventReaderStandard::operator=(const AliFemtoEventReaderStandard& aReader)
 103 {
 104   // Assignment operator
 105   if (this == &aReader)
 106     return *this;
 107
 108   fCurEvent = aReader.fCurEvent;
 109   fCurFile = aReader.fCurFile;
 110   if (fESDEvent) delete fESDEvent;
 111   fESDEvent = new AliESDEvent();
 112   if (fAODEvent) delete fAODEvent;
 113   fAODEvent = new AliAODEvent();
 114   fStack = aReader.fStack;
 115   fGenHeader = aReader.fGenHeader;
 116   fInputType = aReader.fInputType;
 117   fUsePhysicsSel = aReader.fUsePhysicsSel;
 118   if (fUsePhysicsSel) fSelect = new AliPhysicsSelection();
 119   if (fTrackCuts) delete fTrackCuts;
 120   fTrackCuts = new AliESDtrackCuts(*(aReader.fTrackCuts));
 121   fUseTPCOnly = aReader.fUseTPCOnly;
 122
 123   return *this;
 124 }
 125 //__________________
 126 // Simple report
 127 AliFemtoString AliFemtoEventReaderStandard::Report()
 128 {
 129   AliFemtoString temp = "\n This is the AliFemtoEventReaderStandard\n";
 130   return temp;
 131 }
 132 //__________________
 133 AliFemtoEvent* AliFemtoEventReaderStandard::ReturnHbtEvent()
 134 {
 135   // Get the event, read all the relevant information
 136   // and fill the AliFemtoEvent class
 137   // Returns a valid AliFemtoEvent
 138   AliFemtoEvent *hbtEvent = 0;
 139   string tFriendFileName;
 140
 141   hbtEvent = new AliFemtoEvent;
 142
 143   // Get the friend information
 144   cout<<"starting to read event "<<fCurEvent<<endl;
 145   if ((fInputType == kESD) || (fInputType == kESDKine)) {
 146     if(fESDEvent->GetAliESDOld())fESDEvent->CopyFromOldESD();
 147
 148     if (fUsePhysicsSel) {
 149       hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fESDEvent));
 150       if (!(fSelect->IsCollisionCandidate(fESDEvent)))
 151         printf("Event not a collision candidate\n");
 152     }
 153     else
 154       hbtEvent->SetIsCollisionCandidate(kTRUE);
 155   }
 156   else {
 157     hbtEvent->SetIsCollisionCandidate(kTRUE);
 158   }
 159
 160   hbtEvent = new AliFemtoEvent;
 161
 162   double fV1[3];
 163
 164   //setting basic things
 165   if ((fInputType == kESD) || (fInputType == kESDKine)) {
 166     hbtEvent->SetRunNumber(fESDEvent->GetRunNumber());
 167     hbtEvent->SetMagneticField(fESDEvent->GetMagneticField()*kilogauss);//to check if here is ok
 168     hbtEvent->SetZDCN1Energy(fESDEvent->GetZDCN1Energy());
 169     hbtEvent->SetZDCP1Energy(fESDEvent->GetZDCP1Energy());
 170     hbtEvent->SetZDCN2Energy(fESDEvent->GetZDCN2Energy());
 171     hbtEvent->SetZDCP2Energy(fESDEvent->GetZDCP2Energy());
 172     hbtEvent->SetZDCEMEnergy(fESDEvent->GetZDCEMEnergy());
 173     hbtEvent->SetZDCParticipants(fESDEvent->GetZDCParticipants());
 174     hbtEvent->SetTriggerMask(fESDEvent->GetTriggerMask());
 175     hbtEvent->SetTriggerCluster(fESDEvent->GetTriggerCluster());
 176
 177     printf("Got event type %i\n", fESDEvent->GetEventType());
 178
 179     //Vertex
 180     double fVCov[6];
 181     //   if (fUseTPCOnly) {
 182     //     fESDEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 183     //     fESDEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 184     //     if (!fESDEvent->GetPrimaryVertexTPC()->GetStatus())
 185     //       fVCov[4] = -1001.0;
 186     //   }
 187     //   else {
 188     if (fESDEvent->GetPrimaryVertex()) {
 189       fESDEvent->GetPrimaryVertex()->GetXYZ(fV1);
 190       fESDEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 191
 192       if (!fESDEvent->GetPrimaryVertex()->GetStatus()) {
 193         // Get the vertex from SPD
 194         fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
 195         fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
 196
 197
 198         if (!fESDEvent->GetPrimaryVertexSPD()->GetStatus())
 199           fVCov[4] = -1001.0;
 200         else {
 201           fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
 202           fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
 203         }
 204       }
 205     }
 206     else {
 207       if (fESDEvent->GetPrimaryVertexSPD()) {
 208         fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
 209         fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
 210       }
 211     }
 212     if ((!fESDEvent->GetPrimaryVertex()) && (!fESDEvent->GetPrimaryVertexSPD()))
 213       {
 214         cout << "No vertex found !!!" << endl;
 215         fV1[0] = 10000.0;
 216         fV1[1] = 10000.0;
 217         fV1[2] = 10000.0;
 218         fVCov[4] = -1001.0;
 219       }
 220
 221     AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 222
 223     hbtEvent->SetPrimVertPos(vertex);
 224     hbtEvent->SetPrimVertCov(fVCov);
 225   }
 226
 227   if ((fInputType == kAOD) || (fInputType == kAODKine)) {
 228     hbtEvent->SetRunNumber(fAODEvent->GetRunNumber());
 229     hbtEvent->SetMagneticField(fAODEvent->GetMagneticField()*kilogauss);//to check if here is ok
 230     hbtEvent->SetZDCN1Energy(fAODEvent->GetZDCN1Energy());
 231     hbtEvent->SetZDCP1Energy(fAODEvent->GetZDCP1Energy());
 232     hbtEvent->SetZDCN2Energy(fAODEvent->GetZDCN2Energy());
 233     hbtEvent->SetZDCP2Energy(fAODEvent->GetZDCP2Energy());
 234     hbtEvent->SetZDCEMEnergy(fAODEvent->GetZDCEMEnergy(0));
 235     hbtEvent->SetZDCParticipants(0);
 236     hbtEvent->SetTriggerMask(fAODEvent->GetTriggerMask());
 237     hbtEvent->SetTriggerCluster(fAODEvent->GetTriggerCluster());
 238
 239     // Primary Vertex position
 240     fAODEvent->GetPrimaryVertex()->GetPosition(fV1);
 241
 242     AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 243     hbtEvent->SetPrimVertPos(vertex);
 244   }
 245
 246   if ((fInputType == kESDKine) || (fInputType == kAODKine)) {
 247     Double_t tReactionPlane = 0;
 248
 249     AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 250     if (!hdh) {
 251       AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 252       if (cdh) {
 253         TList *tGenHeaders = cdh->GetHeaders();
 254         for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 255           hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 256           if (hdh) break;
 257         }
 258       }
 259     }
 260     if (hdh)
 261       {
 262         tReactionPlane = hdh->ReactionPlaneAngle();
 263         cout << "Got reaction plane " << tReactionPlane << endl;
 264       }
 265
 266     hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 267   }
 268
 269   //starting to reading tracks
 270   int nofTracks=0;  //number of reconstructed tracks in event
 271   if ((fInputType == kESD) || (fInputType == kESDKine))
 272     nofTracks = fESDEvent->GetNumberOfTracks();
 273   else if ((fInputType == kAOD) || (fInputType == kAODKine))
 274     nofTracks = fAODEvent->GetNumberOfTracks();
 275
 276   int realnofTracks=0;//number of track which we use ina analysis
 277
 278   TClonesArray *mcP = 0;
 279   Int_t *motherids=0;
 280
 281   if (fInputType == kAODKine) {
 282     // Attempt to access MC header
 283     AliAODMCHeader *mcH;
 284     mcH = (AliAODMCHeader *) fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
 285     if (!mcH) {
 286       cout << "AOD MC information requested, but no header found!" << endl;
 287     }
 288
 289     mcP = (TClonesArray *) fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
 290     if (!mcP) {
 291       cout << "AOD MC information requested, but no particle array found!" << endl;
 292     }
 293
 294     hbtEvent->SetReactionPlaneAngle(fAODEvent->GetHeader()->GetQTheta(0)/2.0);
 295
 296     if (mcP) {
 297       motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
 298       for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
 299
 300       // Read in mother ids
 301       AliAODMCParticle *motherpart;
 302       for (int ip=0; ip<mcP->GetEntries(); ip++) {
 303         motherpart = (AliAODMCParticle *) mcP->At(ip);
 304         if (motherpart->GetDaughter(0) > 0)
 305           motherids[motherpart->GetDaughter(0)] = ip;
 306         if (motherpart->GetDaughter(1) > 0)
 307           motherids[motherpart->GetDaughter(1)] = ip;
 308       }
 309     }
 310   }
 311
 312   if (fInputType == kESDKine) {
 313     motherids = new Int_t[fStack->GetNtrack()];
 314     for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
 315
 316     // Read in mother ids
 317     TParticle *motherpart;
 318     for (int ip=0; ip<fStack->GetNtrack(); ip++) {
 319       motherpart = fStack->Particle(ip);
 320       if (motherpart->GetDaughter(0) > 0)
 321         motherids[motherpart->GetDaughter(0)] = ip;
 322       if (motherpart->GetDaughter(1) > 0)
 323         motherids[motherpart->GetDaughter(1)] = ip;
 324     }
 325   }
 326
 327   for (int i=0;i<nofTracks;i++)
 328     {
 329       //      cout << "Reading track " << i << endl;
 330       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 331
 332       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 333
 334       if ((fInputType == kESD) || (fInputType == kESDKine)) {
 335
 336         AliESDtrack *esdtrack = 0x0;
 337         if (fUseTPCOnly) {
 338           AliESDtrack *mcp = fESDEvent->GetTrack(i);
 339           esdtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESDEvent, mcp->GetID());
 340           //      printf("Got %p for track %i | ", esdtrack, mcp->GetID());
 341         }
 342         else {
 343           esdtrack = fESDEvent->GetTrack(i);//getting next track
 344         }
 345
 346         if (esdtrack && (fTrackCuts->AcceptTrack(esdtrack))) {
 347
 348           trackCopy->SetCharge((short)esdtrack->GetSign());
 349
 350           //in aliroot we have AliPID
 351           //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 352           //we use only 5 first
 353           double esdpid[5];
 354           esdtrack->GetESDpid(esdpid);
 355           trackCopy->SetPidProbElectron(esdpid[0]);
 356           trackCopy->SetPidProbMuon(esdpid[1]);
 357           trackCopy->SetPidProbPion(esdpid[2]);
 358           trackCopy->SetPidProbKaon(esdpid[3]);
 359           trackCopy->SetPidProbProton(esdpid[4]);
 360
 361           double pxyz[3];
 362           double impact[2];
 363           double covimpact[3];
 364
 365           //      if (fUseTPCOnly) {
 366           //        if (!esdtrack->GetTPCInnerParam()) {
 367           //          cout << "No TPC inner param !" << endl;
 368           //          delete trackCopy;
 369           //          continue;
 370           //        }
 371
 372           //        AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 373           //        param->GetXYZ(rxyz);
 374           //        param->PropagateToDCA(fESDEvent->GetPrimaryVertexTPC(), (fESDEvent->GetMagneticField()), 10000, impact, covimpact);
 375           //        param->GetPxPyPz(pxyz);//reading noconstarined momentum
 376
 377           //        if (fRotateToEventPlane) {
 378           //          double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 379           //          double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 380
 381           //          pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 382           //          pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 383           //        }
 384
 385           //        AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 386           //        if (v.mag() < 0.0001) {
 387           //          //          cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 388           //          delete trackCopy;
 389           //          continue;
 390           //        }
 391
 392           //        trackCopy->SetP(v);//setting momentum
 393           //        trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 394
 395           //        const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 396           //        const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 397           //        //setting helix I do not if it is ok
 398           //        AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fESDEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 399           //        trackCopy->SetHelix(helix);
 400
 401           //        //some stuff which could be useful
 402           //        trackCopy->SetImpactD(impact[0]);
 403           //        trackCopy->SetImpactZ(impact[1]);
 404           //        trackCopy->SetCdd(covimpact[0]);
 405           //        trackCopy->SetCdz(covimpact[1]);
 406           //        trackCopy->SetCzz(covimpact[2]);
 407           //        trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 408
 409           //        delete param;
 410           //      }
 411           //      else {
 412           if (fUseTPCOnly)
 413             tGoodMomentum=esdtrack->GetPxPyPz(pxyz);
 414           else
 415             tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 416           //      printf("Got good momentum %i\n", tGoodMomentum);
 417
 418           AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 419           if (v.Mag() < 0.0001) {
 420
 421             delete trackCopy;
 422             continue;
 423           }
 424
 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)(fESDEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 431           trackCopy->SetHelix(helix);
 432
 433           //some stuff which could be useful
 434           float imp[2];
 435           float cim[3];
 436           esdtrack->GetImpactParameters(imp,cim);
 437
 438           impact[0] = imp[0];
 439           impact[1] = imp[1];
 440           covimpact[0] = cim[0];
 441           covimpact[1] = cim[1];
 442           covimpact[2] = cim[2];
 443
 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(AliESDtrackCuts::GetSigmaToVertex(esdtrack));
 450
 451           trackCopy->SetTrackId(esdtrack->GetID());
 452           trackCopy->SetFlags(esdtrack->GetStatus());
 453           trackCopy->SetLabel(esdtrack->GetLabel());
 454
 455           trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 456           trackCopy->SetITSncls(esdtrack->GetNcls(0));
 457           trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 458           trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 459           trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 460           trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 461           trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 462
 463           trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 464           trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 465
 466           double xtpc[3];
 467           esdtrack->GetInnerXYZ(xtpc);
 468           xtpc[2] -= fV1[2];
 469           trackCopy->SetNominalTPCEntrancePoint(xtpc);
 470
 471           esdtrack->GetOuterXYZ(xtpc);
 472           xtpc[2] -= fV1[2];
 473           trackCopy->SetNominalTPCExitPoint(xtpc);
 474
 475           int indexes[3];
 476           for (int ik=0; ik<3; ik++) {
 477             indexes[ik] = esdtrack->GetKinkIndex(ik);
 478           }
 479           trackCopy->SetKinkIndexes(indexes);
 480
 481           // Fill the hidden information with the simulated data
 482           if (fInputType == kESDKine) {
 483             if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
 484               TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 485
 486               // Check the mother information
 487
 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
 493               //   and is tracked
 494
 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];
 500               fpt = tPart->T();
 501
 502               AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 503               tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 504
 505               fpx *= 1e13;
 506               fpy *= 1e13;
 507               fpz *= 1e13;
 508               fpt *= 1e13;
 509
 510               if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 511                 TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 512                 // Check if this is the same particle stored twice on the stack
 513                 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 514                   // It is the same particle
 515                   // Read in the original freeze-out information
 516                   // and convert it from to [fm]
 517
 518                   // EPOS style
 519                   fpx = mother->Vx()*1e13*0.197327;
 520                   fpy = mother->Vy()*1e13*0.197327;
 521                   fpz = mother->Vz()*1e13*0.197327;
 522                   fpt = mother->T() *1e13*0.197327;
 523
 524
 525                   // Therminator style
 526                   //        fpx = mother->Vx()*1e13;
 527                   //        fpy = mother->Vy()*1e13;
 528                   //        fpz = mother->Vz()*1e13;
 529                   //        fpt = mother->T() *1e13*3e10;
 530
 531                 }
 532               }
 533
 534               tInfo->SetPDGPid(tPart->GetPdgCode());
 535
 536               tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 537               Double_t mass2 = (tPart->Energy() *tPart->Energy() -
 538                                 tPart->Px()*tPart->Px() -
 539                                 tPart->Py()*tPart->Py() -
 540                                 tPart->Pz()*tPart->Pz());
 541               if (mass2>0.0)
 542                 tInfo->SetMass(TMath::Sqrt(mass2));
 543               else
 544                 tInfo->SetMass(0.0);
 545
 546               tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 547               trackCopy->SetHiddenInfo(tInfo);
 548             }
 549             else {
 550               AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 551               tInfo->SetMass(0.0);
 552               double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 553               fpx = fV1[0]*1e13;
 554               fpy = fV1[1]*1e13;
 555               fpz = fV1[2]*1e13;
 556               fpt = 0.0;
 557               tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 558
 559               tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
 560
 561               trackCopy->SetHiddenInfo(tInfo);
 562             }
 563           }
 564           //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl;
 565         }
 566         else
 567           tGoodMomentum = false;
 568
 569         if (fUseTPCOnly)
 570           if (esdtrack) delete esdtrack;
 571       }
 572
 573       if ((fInputType == kAOD) || (fInputType == kAODKine)) {
 574         // Read in the normal AliAODTracks
 575         const AliAODTrack *aodtrack=fAODEvent->GetTrack(i); // getting the AODtrack directly
 576
 577         //      if (!aodtrack->TestFilterBit(fFilterBit))
 578         //        continue;
 579
 580         CopyAODtoFemtoTrack(aodtrack, trackCopy);
 581
 582         if (mcP) {
 583           // Fill the hidden information with the simulated data
 584           //      Int_t pLabel = aodtrack->GetLabel();
 585           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 586
 587           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 588           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 589           if (!tPart) {
 590             fpx = fV1[0];
 591             fpy = fV1[1];
 592             fpz = fV1[2];
 593             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 594             tInfo->SetPDGPid(0);
 595             tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
 596             tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
 597             tInfo->SetMass(0);
 598           }
 599           else {
 600             // Check the mother information
 601
 602             // Using the new way of storing the freeze-out information
 603             // Final state particle is stored twice on the stack
 604             // one copy (mother) is stored with original freeze-out information
 605             //   and is not tracked
 606             // the other one (daughter) is stored with primary vertex position
 607             //   and is tracked
 608
 609             // Freeze-out coordinates
 610             fpx = tPart->Xv() - fV1[0];
 611             fpy = tPart->Yv() - fV1[1];
 612             fpz = tPart->Zv() - fV1[2];
 613             //    fpt = tPart->T();
 614
 615             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 616
 617             fpx *= 1e13;
 618             fpy *= 1e13;
 619             fpz *= 1e13;
 620             //    fpt *= 1e13;
 621
 622             //      cout << "Looking for mother ids " << endl;
 623             if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 624               //        cout << "Got mother id" << endl;
 625               AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 626               // Check if this is the same particle stored twice on the stack
 627               if (mother) {
 628                 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 629                   // It is the same particle
 630                   // Read in the original freeze-out information
 631                   // and convert it from to [fm]
 632
 633                   // EPOS style
 634                   //      fpx = mother->Xv()*1e13*0.197327;
 635                   //      fpy = mother->Yv()*1e13*0.197327;
 636                   //      fpz = mother->Zv()*1e13*0.197327;
 637                   //      fpt = mother->T() *1e13*0.197327*0.5;
 638
 639
 640                   // Therminator style
 641                   fpx = mother->Xv()*1e13;
 642                   fpy = mother->Yv()*1e13;
 643                   fpz = mother->Zv()*1e13;
 644                   //          fpt = mother->T() *1e13*3e10;
 645
 646                 }
 647               }
 648             }
 649
 650             //       if (fRotateToEventPlane) {
 651             //  double tPhi = TMath::ATan2(fpy, fpx);
 652             //  double tRad = TMath::Hypot(fpx, fpy);
 653
 654             //  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 655             //  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 656             //       }
 657
 658             tInfo->SetPDGPid(tPart->GetPdgCode());
 659
 660             //    if (fRotateToEventPlane) {
 661             //      double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 662             //      double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 663
 664             //      tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 665             //                             tRad*TMath::Sin(tPhi - tReactionPlane),
 666             //                             tPart->Pz());
 667             //    }
 668             //       else
 669             tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 670             Double_t mass2 = (tPart->E() *tPart->E() -
 671                               tPart->Px()*tPart->Px() -
 672                               tPart->Py()*tPart->Py() -
 673                               tPart->Pz()*tPart->Pz());
 674             if (mass2>0.0)
 675               tInfo->SetMass(TMath::Sqrt(mass2));
 676             else
 677               tInfo->SetMass(0.0);
 678
 679             tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 680           }
 681           trackCopy->SetHiddenInfo(tInfo);
 682         }
 683
 684         double pxyz[3];
 685         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 686         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 687         // Check the sanity of the tracks - reject zero momentum tracks
 688         if (ktP.Mag() == 0) {
 689           delete trackCopy;
 690           continue;
 691         }
 692
 693
 694       }
 695
 696
 697       //decision if we want this track
 698       //if we using diffrent labels we want that this label was use for first time
 699       //if we use hidden info we want to have match between sim data and ESD
 700       if (tGoodMomentum==true)
 701         {
 702           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 703           realnofTracks++;//real number of tracks
 704           //      delete trackCopy;
 705         }
 706       else
 707         {
 708           delete  trackCopy;
 709         }
 710
 711     }
 712
 713   if (motherids)
 714     delete [] motherids;
 715
 716   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 717   fCurEvent++;
 718   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 719   return hbtEvent;
 720 }
 721 //___________________
 722 void AliFemtoEventReaderStandard::SetESDSource(AliESDEvent *aESD)
 723 {
 724   // The chain loads the ESD for us
 725   // You must provide the address where it can be found
 726   fESDEvent = aESD;
 727 }
 728 //___________________
 729 void AliFemtoEventReaderStandard::SetAODSource(AliAODEvent *aAOD)
 730 {
 731   // The chain loads the ESD for us
 732   // You must provide the address where it can be found
 733   fAODEvent = aAOD;
 734 }
 735 //___________________
 736 void AliFemtoEventReaderStandard::SetStackSource(AliStack *aStack)
 737 {
 738   // The chain loads the stack for us
 739   // You must provide the address where it can be found
 740   fStack = aStack;
 741 }
 742 void AliFemtoEventReaderStandard::SetInputType(AliFemtoInputType aInput)
 743 {
 744   // Set the proper input type
 745   fInputType = aInput;
 746 }
 747 //___________________
 748 void AliFemtoEventReaderStandard::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 749 {
 750   // The chain loads the generator event header for us
 751   // You must provide the address where it can be found
 752   fGenHeader = aGenHeader;
 753 }
 754
 755 void AliFemtoEventReaderStandard::SetUsePhysicsSelection(const bool usephysics)
 756 {
 757   fUsePhysicsSel = usephysics;
 758   if (!fSelect) fSelect = new AliPhysicsSelection();
 759 }
 760
 761 void AliFemtoEventReaderStandard::SetESDTrackCuts(AliESDtrackCuts *esdcuts)
 762 {
 763   // Set external ESD track cuts
 764   fTrackCuts = esdcuts;
 765 }
 766
 767 void AliFemtoEventReaderStandard::SetUseTPCOnly(const bool usetpconly)
 768 {
 769   // Set flag to use TPC only tracks
 770   fUseTPCOnly = usetpconly;
 771 }
 772
 773 void AliFemtoEventReaderStandard::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack,
 774                                                       AliFemtoTrack *tFemtoTrack)
 775 {
 776   // Copy the track information from the AOD into the internal AliFemtoTrack
 777   // If it exists, use the additional information from the PWG2 AOD
 778
 779   // Primary Vertex position
 780   double fV1[3];
 781   fAODEvent->GetPrimaryVertex()->GetPosition(fV1);
 782
 783   tFemtoTrack->SetCharge(tAodTrack->Charge());
 784
 785   //in aliroot we have AliPID
 786   //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 787   //we use only 5 first
 788
 789   // AOD pid has 10 components
 790   double aodpid[10];
 791   tAodTrack->GetPID(aodpid);
 792   tFemtoTrack->SetPidProbElectron(aodpid[0]);
 793   tFemtoTrack->SetPidProbMuon(aodpid[1]);
 794   tFemtoTrack->SetPidProbPion(aodpid[2]);
 795   tFemtoTrack->SetPidProbKaon(aodpid[3]);
 796   tFemtoTrack->SetPidProbProton(aodpid[4]);
 797
 798   double pxyz[3];
 799   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
 800   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 801   tFemtoTrack->SetP(v);//setting momentum
 802   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 803   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 804   //setting track helix
 805   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 806   AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fAODEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
 807   tFemtoTrack->SetHelix(helix);
 808
 809   // Flags
 810   tFemtoTrack->SetTrackId(tAodTrack->GetID());
 811   tFemtoTrack->SetFlags(1);
 812   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
 813
 814   // Track quality information
 815   float covmat[6];
 816   tAodTrack->GetCovMatrix(covmat);
 817   tFemtoTrack->SetImpactD(covmat[0]);
 818   tFemtoTrack->SetImpactZ(covmat[2]);
 819   tFemtoTrack->SetCdd(covmat[3]);
 820   tFemtoTrack->SetCdz(covmat[4]);
 821   tFemtoTrack->SetCzz(covmat[5]);
 822   // This information is only available in the ESD
 823   // We put in fake values or reasonable estimates
 824   tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
 825   tFemtoTrack->SetITSncls(1);
 826   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 827   tFemtoTrack->SetTPCncls(1);
 828   tFemtoTrack->SetTPCnclsF(1);
 829   tFemtoTrack->SetTPCsignalN(1);
 830   tFemtoTrack->SetTPCsignalS(1);
 831
 832   TBits tAllTrue;
 833   TBits tAllFalse;
 834   tAllTrue.ResetAllBits(kTRUE);
 835   tAllFalse.ResetAllBits(kFALSE);
 836
 837
 838   // If not use dummy values
 839   tFemtoTrack->SetTPCClusterMap(tAllTrue);
 840   tFemtoTrack->SetTPCSharedMap(tAllFalse);
 841
 842   double xtpc[3] = {0,0,0};
 843   tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 844   tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 845
 846   int indexes[3];
 847   for (int ik=0; ik<3; ik++) {
 848     indexes[ik] = 0;
 849   }
 850   tFemtoTrack->SetKinkIndexes(indexes);
 851 }
 852
 853 AliAODMCParticle* AliFemtoEventReaderStandard::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
 854 {
 855   if (aLabel < 0) return 0;
 856   AliAODMCParticle *aodP;
 857   Int_t posstack = 0;
 858   if (aLabel > mcP->GetEntries())
 859     posstack = mcP->GetEntries();
 860   else
 861     posstack = aLabel;
 862
 863   aodP = (AliAODMCParticle *) mcP->At(posstack);
 864   if (aodP->GetLabel() > posstack) {
 865     do {
 866       aodP = (AliAODMCParticle *) mcP->At(posstack);
 867       if (aodP->GetLabel() == aLabel) return aodP;
 868       posstack--;
 869     }
 870     while (posstack > 0);
 871   }
 872   else {
 873     do {
 874       aodP = (AliAODMCParticle *) mcP->At(posstack);
 875       if (aodP->GetLabel() == aLabel) return aodP;
 876       posstack++;
 877     }
 878     while (posstack < mcP->GetEntries());
 879   }
 880
 881   return 0;
 882 }