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