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