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
 199
 200
 201       //setting multiplicity
 202         realnofTracks++;//real number of tracks (only primary particles)
 203
 204       //setting normalized multiplicity
 205
 206         if(kinetrack->GetPDG()->Charge()/3!=0)
 207           if (kinetrack->Pt() > 0.15 && kinetrack->Pt() < 20)
 208             if (kinetrack->Eta() < 0.8)
 209               tNormMult++;
 210
 211         //counting particles that go into direction of VZERO detector
 212         if(kinetrack->Eta() > 2.8 && kinetrack->Eta() < 5.1) //VZERO-A
 213           tV0direction++;
 214         if(kinetrack->Eta() > -3.7 && kinetrack->Eta() < -1.7)//VZERO-C
 215           tV0direction++;
 216
 217           //charge
 218       trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
 219
 220
 221       //in aliroot we have AliPID
 222       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 223       //we use only 5 first
 224       double kinepid[5];
 225       for(int pid_iter=0;pid_iter<5;pid_iter++)
 226           kinepid[pid_iter]=0;
 227
 228       int pdgcode = kinetrack->GetPdgCode();
 229       //proton
 230       if(pdgcode==2212 || pdgcode==-2212)
 231         kinepid[4]=1000;
 232       //kaon
 233       else if(pdgcode==321 || pdgcode==-321 )
 234         kinepid[3]=1000;
 235       //pion
 236       else if( pdgcode==211 || pdgcode==-211)
 237         kinepid[2]=1000;
 238       //electron
 239       else if(pdgcode==11 || pdgcode==-11)
 240         kinepid[0]=1000;
 241       //muon
 242       else if(pdgcode==13 || pdgcode==-13)
 243         kinepid[1]=1000;
 244       else
 245         continue;
 246       trackCopy->SetPidProbElectron(kinepid[0]);
 247       trackCopy->SetPidProbMuon(kinepid[1]);
 248       trackCopy->SetPidProbPion(kinepid[2]);
 249       trackCopy->SetPidProbKaon(kinepid[3]);
 250       trackCopy->SetPidProbProton(kinepid[4]);
 251
 252
 253         //Momentum
 254       double pxyz[3];
 255       double rxyz[3];
 256
 257       pxyz[0]=kinetrack->Px();
 258       pxyz[1]=kinetrack->Py();
 259       pxyz[2]=kinetrack->Pz();
 260
 261       rxyz[0]=kinetrack->Vx();
 262       rxyz[1]=kinetrack->Vy();
 263       rxyz[2]=kinetrack->Vz();
 264
 265       AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 266       tInfo->SetPDGPid(kinetrack->GetPdgCode());
 267       tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 268       tInfo->SetMass(kinetrack->GetMass());
 269       tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 270       trackCopy->SetHiddenInfo(tInfo);
 271
 272       trackCopy->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 273       trackCopy->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 274
 275
 276         if (fRotateToEventPlane) {
 277           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 278           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 279
 280           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 281           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 282         }
 283
 284         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 285         if (v.Mag() < 0.0001) {
 286           //cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 287           delete trackCopy;
 288           continue;
 289         }
 290
 291         trackCopy->SetP(v);//setting momentum
 292         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 293         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 294         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 295
 296         //label
 297         trackCopy->SetLabel(i);
 298
 299
 300         hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 301         //cout<<"Track added: "<<i<<endl;
 302
 303     }
 304
 305   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 306   if (fEstEventMult == kGlobalCount)
 307     hbtEvent->SetNormalizedMult(tNormMult);
 308   else if(fEstEventMult == kVZERO)
 309     hbtEvent->SetNormalizedMult(tV0direction);
 310   fCurEvent++;
 311
 312
 313   //V0 analysis code - no V0 finder for Kinematics, we can only check if it is primary and if it has at least 2 daughters.
 314
 315   for (int i=0;i<nofTracks;i++)
 316     {
 317       //do not take primaries
 318       if(!fStack->IsPhysicalPrimary(i)) {continue;}
 319       //getting next track
 320       TParticle *kinetrack= fStack->Particle(i);
 321       if (!kinetrack) continue;
 322
 323       if(kinetrack->GetPDG()->Charge()!=0) continue; //charge - neutral
 324       //if(kinetrack->GetPDG()->Stable()==1) continue; //particle is not stable
 325       if(kinetrack->GetDaughter(0)<1) continue; //has 1'st daughter
 326       if(kinetrack->GetDaughter(1)<1) continue;  //has 2'nd daughter
 327
 328       //we want one positive, one negative particle. Or two neutral.
 329       // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()>=0)
 330       //        if((fStack->Particle(kinetrack->GetDaughter(1)))->GetPDG()->Charge()>0)
 331       //          continue;
 332       // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<=0)
 333       //        if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<0)
 334       //          continue;
 335
 336       if(kinetrack->Pt()<0.00001)
 337         continue;
 338
 339       AliFemtoV0* trackCopyV0 = new AliFemtoV0();
 340       CopyAODtoFemtoV0(kinetrack, trackCopyV0);
 341       hbtEvent->V0Collection()->push_back(trackCopyV0);
 342     //cout<<"Pushback v0 to v0collection"<<endl;
 343     }
 344
 345
 346   cout<<"Number of tracks: "<<realnofTracks<<endl;
 347
 348   return hbtEvent;
 349 }
 350
 351 //___________________
 352 void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
 353 {
 354   // The chain loads the stack for us
 355   // You must provide the address where it can be found
 356   fStack = aStack;
 357 }
 358 //___________________
 359 void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 360 {
 361   // The chain loads the generator event header for us
 362   // You must provide the address where it can be found
 363   fGenHeader = aGenHeader;
 364 }
 365
 366 //__________________
 367 void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
 368 {
 369   fRotateToEventPlane=dorotate;
 370 }
 371
 372 void AliFemtoEventReaderKinematicsChain::SetUseMultiplicity(EstEventMult aType)
 373 {
 374   fEstEventMult = aType;
 375 }
 376
 377 Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double *impact, double *covar)
 378 {
 379   // Calculates the number of sigma to the vertex.
 380
 381   Float_t b[2];
 382   Float_t bRes[2];
 383   Float_t bCov[3];
 384
 385   b[0] = impact[0];
 386   b[1] = impact[1];
 387   bCov[0] = covar[0];
 388   bCov[1] = covar[1];
 389   bCov[2] = covar[2];
 390
 391   bRes[0] = TMath::Sqrt(bCov[0]);
 392   bRes[1] = TMath::Sqrt(bCov[2]);
 393
 394   // -----------------------------------
 395   // How to get to a n-sigma cut?
 396   //
 397   // The accumulated statistics from 0 to d is
 398   //
 399   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 400   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 401   //
 402   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 403   // Can this be expressed in a different way?
 404
 405   if (bRes[0] == 0 || bRes[1] ==0)
 406     return -1;
 407
 408   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 409
 410   // stupid rounding problem screws up everything:
 411   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 412   if (TMath::Exp(-d * d / 2) < 1e-10)
 413     return 1000;
 414
 415   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 416   return d;
 417 }
 418
 419
 420
 421  void AliFemtoEventReaderKinematicsChain::CopyAODtoFemtoV0(TParticle *tv0, AliFemtoV0 *tFemtoV0 )
 422 {
 423   tFemtoV0->SetEtaV0(tv0->Eta());
 424   tFemtoV0->SetEtaV0(tv0->Phi());
 425   tFemtoV0->SetptV0(tv0->Pt());
 426   tFemtoV0->SetptotV0(tv0->P());
 427
 428   tFemtoV0->SetmomV0X(tv0->Px());
 429   tFemtoV0->SetmomV0Y(tv0->Py());
 430   tFemtoV0->SetmomV0Z(tv0->Pz());
 431   AliFemtoThreeVector momv0(tv0->Px(),tv0->Py(),tv0->Pz());
 432   tFemtoV0->SetmomV0(momv0);
 433
 434
 435   TParticle *trackpos;
 436   TParticle *trackneg;
 437
 438   //daughters
 439   if(fStack->Particle(tv0->GetDaughter(0))->GetPDG()->Charge()>=0) //first positive, second negative
 440     {
 441       trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
 442       trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
 443       tFemtoV0->SetidPos(tv0->GetDaughter(0));
 444       tFemtoV0->SetidNeg(tv0->GetDaughter(1));
 445     }
 446   else //first negative, second positive
 447     {
 448       trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
 449       trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
 450       tFemtoV0->SetidPos(tv0->GetDaughter(1));
 451       tFemtoV0->SetidNeg(tv0->GetDaughter(0));
 452     }
 453
 454   tFemtoV0->SetEtaPos(trackpos->Eta());
 455   tFemtoV0->SetEtaNeg(trackneg->Eta());
 456
 457   tFemtoV0->SetptPos(trackpos->Pt());
 458   tFemtoV0->SetptNeg(trackneg->Pt());
 459
 460   tFemtoV0->SetptotPos(trackpos->P());
 461   tFemtoV0->SetptotNeg(trackneg->P());
 462
 463   tFemtoV0->SetmomPosX(trackpos->Px());
 464   tFemtoV0->SetmomPosY(trackpos->Py());
 465   tFemtoV0->SetmomPosZ(trackpos->Pz());
 466   AliFemtoThreeVector mompos(trackpos->Px(),trackpos->Py(),trackpos->Pz());
 467   tFemtoV0->SetmomPos(mompos);
 468
 469   tFemtoV0->SetmomNegX(trackneg->Px());
 470   tFemtoV0->SetmomNegY(trackneg->Py());
 471   tFemtoV0->SetmomNegZ(trackneg->Pz());
 472   AliFemtoThreeVector momneg(trackneg->Px(),trackneg->Py(),trackneg->Pz());
 473   tFemtoV0->SetmomNeg(momneg);
 474
 475
 476   tFemtoV0->SetmassLambda(tv0->GetMass());
 477   tFemtoV0->SetmassAntiLambda(tv0->GetMass());
 478   tFemtoV0->SetmassK0Short(tv0->GetMass());
 479
 480   tFemtoV0->SetYV0(tv0->Y());
 481
 482   tFemtoV0->SetdecayVertexV0X(trackpos->Vx()); //vertex of the decay is set as the vertex of creation of daughters
 483   tFemtoV0->SetdecayVertexV0Y(trackpos->Vy());
 484   tFemtoV0->SetdecayVertexV0Z(trackpos->Vz());
 485   AliFemtoThreeVector decayvertex(trackpos->Vx(),trackpos->Vy(),trackpos->Vz());
 486   tFemtoV0->SetdecayVertexV0(decayvertex);
 487
 488   tFemtoV0->SetdcaV0Daughters(0);
 489   tFemtoV0->SetCosPointingAngle(1);
 490
 491
 492   tFemtoV0->SetStatusPos(1);
 493   tFemtoV0->SetStatusNeg(1);
 494
 495
 496   if(trackpos->GetPdgCode()==2212) //proton
 497     {
 498       tFemtoV0->SetPosNSigmaTPCK(1000);
 499       tFemtoV0->SetPosNSigmaTPCPi(1000);
 500       tFemtoV0->SetPosNSigmaTPCP(0);
 501     }
 502   if(trackneg->GetPdgCode()==-2212) //antiproton
 503     {
 504       tFemtoV0->SetNegNSigmaTPCK(1000);
 505       tFemtoV0->SetNegNSigmaTPCPi(1000);
 506       tFemtoV0->SetNegNSigmaTPCP(0);
 507     }
 508   if(trackpos->GetPdgCode()==211) //pion plus
 509     {
 510       tFemtoV0->SetPosNSigmaTPCK(1000);
 511       tFemtoV0->SetPosNSigmaTPCPi(0);
 512       tFemtoV0->SetPosNSigmaTPCP(1000);
 513     }
 514   if(trackneg->GetPdgCode()==-211) //pion minus
 515     {
 516       tFemtoV0->SetNegNSigmaTPCK(1000);
 517       tFemtoV0->SetNegNSigmaTPCPi(0);
 518       tFemtoV0->SetNegNSigmaTPCP(1000);
 519     }
 520   if(trackpos->GetPdgCode()==321) //K+
 521     {
 522       tFemtoV0->SetPosNSigmaTPCK(0);
 523       tFemtoV0->SetPosNSigmaTPCPi(1000);
 524       tFemtoV0->SetPosNSigmaTPCP(1000);
 525     }
 526   if(trackneg->GetPdgCode()==-321) //K-
 527     {
 528       tFemtoV0->SetNegNSigmaTPCK(0);
 529       tFemtoV0->SetNegNSigmaTPCPi(1000);
 530       tFemtoV0->SetNegNSigmaTPCP(1000);
 531     }
 532
 533
 534 }