PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.cxx

   1 /////////////////////////////////////////////////////////////////////////////////////
   2 //                                                                                 //
   3 // AliFemtoEventReaderKinematicsChain - the reader class for the Alice ESD and     //
   4 // the model Kinematics information tailored for the Task framework and the        //
   5 // Reads in AliESDfriend to create shared hit/quality information                  //
   6 // Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
   7 //          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
   8 //                                                                                 //
   9 /////////////////////////////////////////////////////////////////////////////////////
  10
  11 #include "AliFemtoEventReaderKinematicsChain.h"
  12
  13 #include "TFile.h"
  14 #include "TTree.h"
  15 #include "TList.h"
  16
  17 #include "AliFmPhysicalHelixD.h"
  18 #include "AliFmThreeVectorF.h"
  19
  20 #include "SystemOfUnits.h"
  21
  22 #include "AliFemtoEvent.h"
  23
  24 #include "TParticle.h"
  25 #include "AliStack.h"
  26 #include "TParticlePDG.h"
  27 #include "AliFemtoModelHiddenInfo.h"
  28 #include "AliFemtoModelGlobalHiddenInfo.h"
  29 #include "AliGenHijingEventHeader.h"
  30 #include "AliGenCocktailEventHeader.h"
  31
  32 #include "AliVertexerTracks.h"
  33
  34 ClassImp(AliFemtoEventReaderKinematicsChain)
  35
  36 #if !(ST_NO_NAMESPACES)
  37   using namespace units;
  38 #endif
  39
  40 using namespace std;
  41 //____________________________
  42 AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain():
  43   fFileName(" "),
  44   fConstrained(true),
  45   fNumberofEvent(0),
  46   fCurEvent(0),
  47   fCurFile(0),
  48   fStack(0x0),
  49   fGenHeader(0x0),
  50   fEstEventMult(kGlobalCount),
  51   fRotateToEventPlane(0)
  52 {
  53   //constructor with 0 parameters , look at default settings
  54 }
  55
  56 //__________________
  57 AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
  58   AliFemtoEventReader(aReader),
  59   fFileName(" "),
  60   fConstrained(true),
  61   fNumberofEvent(0),
  62   fCurEvent(0),
  63   fCurFile(0),
  64   fStack(0x0),
  65   fGenHeader(0x0),
  66   fEstEventMult(kGlobalCount),
  67   fRotateToEventPlane(0)
  68 {
  69   // Copy constructor
  70   fConstrained = aReader.fConstrained;
  71   fNumberofEvent = aReader.fNumberofEvent;
  72   fCurEvent = aReader.fCurEvent;
  73   fCurFile = aReader.fCurFile;
  74   fStack = aReader.fStack;
  75   fEstEventMult = aReader.fEstEventMult;
  76   fRotateToEventPlane = aReader.fRotateToEventPlane;
  77 }
  78 //__________________
  79 AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
  80 {
  81   //Destructor
  82   //delete fEvent;
  83 }
  84
  85 //__________________
  86 AliFemtoEventReaderKinematicsChain& AliFemtoEventReaderKinematicsChain::operator=(const AliFemtoEventReaderKinematicsChain& aReader)
  87 {
  88   // Assignment operator
  89   if (this == &aReader)
  90     return *this;
  91
  92   fConstrained = aReader.fConstrained;
  93   fNumberofEvent = aReader.fNumberofEvent;
  94   fCurEvent = aReader.fCurEvent;
  95   fCurFile = aReader.fCurFile;
  96   fStack = aReader.fStack;
  97   fGenHeader = aReader.fGenHeader;
  98   fEstEventMult = aReader.fEstEventMult;
  99   fRotateToEventPlane = aReader.fRotateToEventPlane;
 100   return *this;
 101 }
 102 //__________________
 103 // Simple report
 104 AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
 105 {
 106   AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
 107   return temp;
 108 }
 109
 110 //__________________
 111 void AliFemtoEventReaderKinematicsChain::SetConstrained(const bool constrained)
 112 {
 113   // Select whether to read constrained or not constrained momentum
 114   fConstrained=constrained;
 115 }
 116 //__________________
 117 bool AliFemtoEventReaderKinematicsChain::GetConstrained() const
 118 {
 119   // Check whether we read constrained or not constrained momentum
 120   return fConstrained;
 121 }
 122 //__________________
 123 AliFemtoEvent* AliFemtoEventReaderKinematicsChain::ReturnHbtEvent()
 124 {
 125   // Get the event, read all the relevant information from the stack
 126   // and fill the AliFemtoEvent class
 127   // Returns a valid AliFemtoEvent
 128   AliFemtoEvent *hbtEvent = 0;
 129   string tFriendFileName;
 130
 131   cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
 132
 133   hbtEvent = new AliFemtoEvent;
 134   //setting basic things
 135   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 136   hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
 137   hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
 138   hbtEvent->SetZDCN1Energy(0);
 139   hbtEvent->SetZDCP1Energy(0);
 140   hbtEvent->SetZDCN2Energy(0);
 141   hbtEvent->SetZDCP2Energy(0);
 142   hbtEvent->SetZDCEMEnergy(0);
 143   hbtEvent->SetZDCParticipants(0);
 144   hbtEvent->SetTriggerMask(0);
 145   hbtEvent->SetTriggerCluster(0);
 146
 147   //Vertex
 148   double fV1[3] = {0.0,0.0,0.0};
 149   double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 150
 151
 152   AliFmThreeVectorF vertex(0,0,0);
 153
 154
 155   hbtEvent->SetPrimVertPos(vertex);
 156   hbtEvent->SetPrimVertCov(fVCov);
 157
 158   Double_t tReactionPlane = 0;
 159
 160   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 161   if (!hdh) {
 162     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 163     if (cdh) {
 164       TList *tGenHeaders = cdh->GetHeaders();
 165       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 166         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 167         if (hdh) break;
 168       }
 169     }
 170   }
 171
 172   if (hdh)
 173     {
 174       tReactionPlane = hdh->ReactionPlaneAngle();
 175       cout << "Got reaction plane " << tReactionPlane << endl;
 176     }
 177
 178   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 179
 180   //starting to reading tracks
 181   int nofTracks=0;  //number of all tracks in MC event
 182   nofTracks=fStack->GetNtrack();
 183   int realnofTracks=0;//number of track which we use in analysis
 184
 185
 186   int tNormMult = 0;
 187   int tV0direction = 0;
 188   for (int i=0;i<nofTracks;i++)
 189     {
 190       //take only primaries
 191       if(!fStack->IsPhysicalPrimary(i)) {continue;}
 192
 193       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 194
 195           //getting next track
 196       TParticle *kinetrack= fStack->Particle(i);
 197
 198       //setting multiplicity
 199         realnofTracks++;//real number of tracks (only primary particles)
 200
 201       //setting normalized multiplicity
 202
 203         if(kinetrack->GetPDG()->Charge()/3!=0)
 204           if (kinetrack->Pt() > 0.15 && kinetrack->Pt() < 20)
 205             if (kinetrack->Eta() < 0.8)
 206               tNormMult++;
 207
 208         //counting particles that go into direction of VZERO detector
 209         if(kinetrack->Eta() > 2.8 && kinetrack->Eta() < 5.1) //VZERO-A
 210           tV0direction++;
 211         if(kinetrack->Eta() > -3.7 && kinetrack->Eta() < -1.7)//VZERO-C
 212           tV0direction++;
 213
 214           //charge
 215       trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
 216
 217
 218       //in aliroot we have AliPID
 219       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 220       //we use only 5 first
 221       double kinepid[5];
 222       for(int pid_iter=0;pid_iter<5;pid_iter++)
 223           kinepid[pid_iter]=0;
 224
 225       int pdgcode = kinetrack->GetPdgCode();
 226       //proton
 227       if(pdgcode==2212 || pdgcode==-2212)
 228         kinepid[4]=1000;
 229       //kaon
 230       else if(pdgcode==321 || pdgcode==-321 )
 231         kinepid[3]=1000;
 232       //pion
 233       else if( pdgcode==211 || pdgcode==-211)
 234         kinepid[2]=1000;
 235       //electron
 236       else if(pdgcode==11 || pdgcode==-11)
 237         kinepid[0]=1000;
 238       //muon
 239       else if(pdgcode==13 || pdgcode==-13)
 240         kinepid[1]=1000;
 241       else
 242         continue;
 243       trackCopy->SetPidProbElectron(kinepid[0]);
 244       trackCopy->SetPidProbMuon(kinepid[1]);
 245       trackCopy->SetPidProbPion(kinepid[2]);
 246       trackCopy->SetPidProbKaon(kinepid[3]);
 247       trackCopy->SetPidProbProton(kinepid[4]);
 248
 249
 250         //Momentum
 251       double pxyz[3];
 252       // double rxyz[3];
 253
 254         pxyz[0]=kinetrack->Px();
 255         pxyz[1]=kinetrack->Py();
 256         pxyz[2]=kinetrack->Pz();
 257
 258         // rxyz[0]=kinetrack->Vx();
 259         // rxyz[1]=kinetrack->Vy();
 260         // rxyz[2]=kinetrack->Vz();
 261
 262         if (fRotateToEventPlane) {
 263           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 264           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 265
 266           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 267           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 268         }
 269
 270         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 271         if (v.Mag() < 0.0001) {
 272           //cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 273           delete trackCopy;
 274           continue;
 275         }
 276
 277         trackCopy->SetP(v);//setting momentum
 278         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 279         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 280         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 281
 282         //label
 283         trackCopy->SetLabel(i);
 284
 285
 286         hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 287         //cout<<"Track added: "<<i<<endl;
 288
 289     }
 290
 291   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 292   if (fEstEventMult == kGlobalCount)
 293     hbtEvent->SetNormalizedMult(tNormMult);
 294   else if(fEstEventMult == kVZERO)
 295     hbtEvent->SetNormalizedMult(tV0direction);
 296   fCurEvent++;
 297
 298
 299   //V0 analysis code - no V0 finder for Kinematics, we can only check if it is primary and if it has at least 2 daughters.
 300
 301   for (int i=0;i<nofTracks;i++)
 302     {
 303       //do not take primaries
 304       if(!fStack->IsPhysicalPrimary(i)) {continue;}
 305       //getting next track
 306       TParticle *kinetrack= fStack->Particle(i);
 307       if (!kinetrack) continue;
 308
 309       if(kinetrack->GetPDG()->Charge()!=0) continue; //charge - neutral
 310       //if(kinetrack->GetPDG()->Stable()==1) continue; //particle is not stable
 311       if(kinetrack->GetDaughter(0)<1) continue; //has 1'st daughter
 312       if(kinetrack->GetDaughter(1)<1) continue;  //has 2'nd daughter
 313
 314       //we want one positive, one negative particle. Or two neutral.
 315       // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()>=0)
 316       //        if((fStack->Particle(kinetrack->GetDaughter(1)))->GetPDG()->Charge()>0)
 317       //          continue;
 318       // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<=0)
 319       //        if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<0)
 320       //          continue;
 321
 322       if(kinetrack->Pt()<0.00001)
 323         continue;
 324
 325       AliFemtoV0* trackCopyV0 = new AliFemtoV0();
 326       CopyAODtoFemtoV0(kinetrack, trackCopyV0);
 327       hbtEvent->V0Collection()->push_back(trackCopyV0);
 328     //cout<<"Pushback v0 to v0collection"<<endl;
 329     }
 330
 331
 332   cout<<"Number of tracks: "<<realnofTracks<<endl;
 333
 334   return hbtEvent;
 335 }
 336
 337 //___________________
 338 void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
 339 {
 340   // The chain loads the stack for us
 341   // You must provide the address where it can be found
 342   fStack = aStack;
 343 }
 344 //___________________
 345 void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 346 {
 347   // The chain loads the generator event header for us
 348   // You must provide the address where it can be found
 349   fGenHeader = aGenHeader;
 350 }
 351
 352 //__________________
 353 void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
 354 {
 355   fRotateToEventPlane=dorotate;
 356 }
 357
 358 void AliFemtoEventReaderKinematicsChain::SetUseMultiplicity(EstEventMult aType)
 359 {
 360   fEstEventMult = aType;
 361 }
 362
 363 Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double *impact, double *covar)
 364 {
 365   // Calculates the number of sigma to the vertex.
 366
 367   Float_t b[2];
 368   Float_t bRes[2];
 369   Float_t bCov[3];
 370
 371   b[0] = impact[0];
 372   b[1] = impact[1];
 373   bCov[0] = covar[0];
 374   bCov[1] = covar[1];
 375   bCov[2] = covar[2];
 376
 377   bRes[0] = TMath::Sqrt(bCov[0]);
 378   bRes[1] = TMath::Sqrt(bCov[2]);
 379
 380   // -----------------------------------
 381   // How to get to a n-sigma cut?
 382   //
 383   // The accumulated statistics from 0 to d is
 384   //
 385   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 386   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 387   //
 388   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 389   // Can this be expressed in a different way?
 390
 391   if (bRes[0] == 0 || bRes[1] ==0)
 392     return -1;
 393
 394   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 395
 396   // stupid rounding problem screws up everything:
 397   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 398   if (TMath::Exp(-d * d / 2) < 1e-10)
 399     return 1000;
 400
 401   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 402   return d;
 403 }
 404
 405
 406
 407  void AliFemtoEventReaderKinematicsChain::CopyAODtoFemtoV0(TParticle *tv0, AliFemtoV0 *tFemtoV0 )
 408 {
 409   tFemtoV0->SetEtaV0(tv0->Eta());
 410   tFemtoV0->SetEtaV0(tv0->Phi());
 411   tFemtoV0->SetptV0(tv0->Pt());
 412   tFemtoV0->SetptotV0(tv0->P());
 413
 414   tFemtoV0->SetmomV0X(tv0->Px());
 415   tFemtoV0->SetmomV0Y(tv0->Py());
 416   tFemtoV0->SetmomV0Z(tv0->Pz());
 417   AliFemtoThreeVector momv0(tv0->Px(),tv0->Py(),tv0->Pz());
 418   tFemtoV0->SetmomV0(momv0);
 419
 420
 421   TParticle *trackpos;
 422   TParticle *trackneg;
 423
 424   //daughters
 425   if(fStack->Particle(tv0->GetDaughter(0))->GetPDG()->Charge()>=0) //first positive, second negative
 426     {
 427       trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
 428       trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
 429       tFemtoV0->SetidPos(tv0->GetDaughter(0));
 430       tFemtoV0->SetidNeg(tv0->GetDaughter(1));
 431     }
 432   else //first negative, second positive
 433     {
 434       trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
 435       trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
 436       tFemtoV0->SetidPos(tv0->GetDaughter(1));
 437       tFemtoV0->SetidNeg(tv0->GetDaughter(0));
 438     }
 439
 440   tFemtoV0->SetEtaPos(trackpos->Eta());
 441   tFemtoV0->SetEtaNeg(trackneg->Eta());
 442
 443   tFemtoV0->SetptPos(trackpos->Pt());
 444   tFemtoV0->SetptNeg(trackneg->Pt());
 445
 446   tFemtoV0->SetptotPos(trackpos->P());
 447   tFemtoV0->SetptotNeg(trackneg->P());
 448
 449   tFemtoV0->SetmomPosX(trackpos->Px());
 450   tFemtoV0->SetmomPosY(trackpos->Py());
 451   tFemtoV0->SetmomPosZ(trackpos->Pz());
 452   AliFemtoThreeVector mompos(trackpos->Px(),trackpos->Py(),trackpos->Pz());
 453   tFemtoV0->SetmomPos(mompos);
 454
 455   tFemtoV0->SetmomNegX(trackneg->Px());
 456   tFemtoV0->SetmomNegY(trackneg->Py());
 457   tFemtoV0->SetmomNegZ(trackneg->Pz());
 458   AliFemtoThreeVector momneg(trackneg->Px(),trackneg->Py(),trackneg->Pz());
 459   tFemtoV0->SetmomNeg(momneg);
 460
 461
 462   tFemtoV0->SetmassLambda(tv0->GetMass());
 463   tFemtoV0->SetmassAntiLambda(tv0->GetMass());
 464   tFemtoV0->SetmassK0Short(tv0->GetMass());
 465
 466   tFemtoV0->SetYV0(tv0->Y());
 467
 468   tFemtoV0->SetdecayVertexV0X(trackpos->Vx()); //vertex of the decay is set as the vertex of creation of daughters
 469   tFemtoV0->SetdecayVertexV0Y(trackpos->Vy());
 470   tFemtoV0->SetdecayVertexV0Z(trackpos->Vz());
 471   AliFemtoThreeVector decayvertex(trackpos->Vx(),trackpos->Vy(),trackpos->Vz());
 472   tFemtoV0->SetdecayVertexV0(decayvertex);
 473
 474   tFemtoV0->SetdcaV0Daughters(0);
 475   tFemtoV0->SetCosPointingAngle(1);
 476
 477
 478   tFemtoV0->SetStatusPos(1);
 479   tFemtoV0->SetStatusNeg(1);
 480
 481
 482   if(trackpos->GetPdgCode()==2212) //proton
 483     {
 484       tFemtoV0->SetPosNSigmaTPCK(1000);
 485       tFemtoV0->SetPosNSigmaTPCPi(1000);
 486       tFemtoV0->SetPosNSigmaTPCP(0);
 487     }
 488   if(trackneg->GetPdgCode()==-2212) //antiproton
 489     {
 490       tFemtoV0->SetNegNSigmaTPCK(1000);
 491       tFemtoV0->SetNegNSigmaTPCPi(1000);
 492       tFemtoV0->SetNegNSigmaTPCP(0);
 493     }
 494   if(trackpos->GetPdgCode()==211) //pion plus
 495     {
 496       tFemtoV0->SetPosNSigmaTPCK(1000);
 497       tFemtoV0->SetPosNSigmaTPCPi(0);
 498       tFemtoV0->SetPosNSigmaTPCP(1000);
 499     }
 500   if(trackneg->GetPdgCode()==-211) //pion minus
 501     {
 502       tFemtoV0->SetNegNSigmaTPCK(1000);
 503       tFemtoV0->SetNegNSigmaTPCPi(0);
 504       tFemtoV0->SetNegNSigmaTPCP(1000);
 505     }
 506   if(trackpos->GetPdgCode()==321) //K+
 507     {
 508       tFemtoV0->SetPosNSigmaTPCK(0);
 509       tFemtoV0->SetPosNSigmaTPCPi(1000);
 510       tFemtoV0->SetPosNSigmaTPCP(1000);
 511     }
 512   if(trackneg->GetPdgCode()==-321) //K-
 513     {
 514       tFemtoV0->SetNegNSigmaTPCK(0);
 515       tFemtoV0->SetNegNSigmaTPCPi(1000);
 516       tFemtoV0->SetNegNSigmaTPCP(1000);
 517     }
 518
 519
 520 }