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   fRotateToEventPlane(0)
  51 {
  52   //constructor with 0 parameters , look at default settings
  53 }
  54
  55 //__________________
  56 AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
  57   AliFemtoEventReader(aReader),
  58   fFileName(" "),
  59   fConstrained(true),
  60   fNumberofEvent(0),
  61   fCurEvent(0),
  62   fCurFile(0),
  63   fStack(0x0),
  64   fGenHeader(0x0),
  65   fRotateToEventPlane(0)
  66 {
  67   // Copy constructor
  68   fConstrained = aReader.fConstrained;
  69   fNumberofEvent = aReader.fNumberofEvent;
  70   fCurEvent = aReader.fCurEvent;
  71   fCurFile = aReader.fCurFile;
  72   fStack = aReader.fStack;
  73   fRotateToEventPlane = aReader.fRotateToEventPlane;
  74 }
  75 //__________________
  76 AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
  77 {
  78   //Destructor
  79   //delete fEvent;
  80 }
  81
  82 //__________________
  83 AliFemtoEventReaderKinematicsChain& AliFemtoEventReaderKinematicsChain::operator=(const AliFemtoEventReaderKinematicsChain& aReader)
  84 {
  85   // Assignment operator
  86   if (this == &aReader)
  87     return *this;
  88
  89   fConstrained = aReader.fConstrained;
  90   fNumberofEvent = aReader.fNumberofEvent;
  91   fCurEvent = aReader.fCurEvent;
  92   fCurFile = aReader.fCurFile;
  93   fStack = aReader.fStack;
  94   fGenHeader = aReader.fGenHeader;
  95   fRotateToEventPlane = aReader.fRotateToEventPlane;
  96   return *this;
  97 }
  98 //__________________
  99 // Simple report
 100 AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
 101 {
 102   AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
 103   return temp;
 104 }
 105
 106 //__________________
 107 void AliFemtoEventReaderKinematicsChain::SetConstrained(const bool constrained)
 108 {
 109   // Select whether to read constrained or not constrained momentum
 110   fConstrained=constrained;
 111 }
 112 //__________________
 113 bool AliFemtoEventReaderKinematicsChain::GetConstrained() const
 114 {
 115   // Check whether we read constrained or not constrained momentum
 116   return fConstrained;
 117 }
 118 //__________________
 119 AliFemtoEvent* AliFemtoEventReaderKinematicsChain::ReturnHbtEvent()
 120 {
 121   // Get the event, read all the relevant information from the stack
 122   // and fill the AliFemtoEvent class
 123   // Returns a valid AliFemtoEvent
 124   AliFemtoEvent *hbtEvent = 0;
 125   string tFriendFileName;
 126
 127   cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
 128
 129   hbtEvent = new AliFemtoEvent;
 130   //setting basic things
 131   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 132   hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
 133   hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
 134   hbtEvent->SetZDCN1Energy(0);
 135   hbtEvent->SetZDCP1Energy(0);
 136   hbtEvent->SetZDCN2Energy(0);
 137   hbtEvent->SetZDCP2Energy(0);
 138   hbtEvent->SetZDCEMEnergy(0);
 139   hbtEvent->SetZDCParticipants(0);
 140   hbtEvent->SetTriggerMask(0);
 141   hbtEvent->SetTriggerCluster(0);
 142
 143   //Vertex
 144   double fV1[3] = {0.0,0.0,0.0};
 145   double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 146
 147
 148   AliFmThreeVectorF vertex(0,0,0);
 149
 150
 151   hbtEvent->SetPrimVertPos(vertex);
 152   hbtEvent->SetPrimVertCov(fVCov);
 153
 154   Double_t tReactionPlane = 0;
 155
 156   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 157   if (!hdh) {
 158     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 159     if (cdh) {
 160       TList *tGenHeaders = cdh->GetHeaders();
 161       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 162         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 163         if (hdh) break;
 164       }
 165     }
 166   }
 167
 168   if (hdh)
 169     {
 170       tReactionPlane = hdh->ReactionPlaneAngle();
 171       cout << "Got reaction plane " << tReactionPlane << endl;
 172     }
 173
 174   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 175
 176   //starting to reading tracks
 177   int nofTracks=0;  //number of all tracks in MC event
 178   nofTracks=fStack->GetNtrack();
 179   int realnofTracks=0;//number of track which we use in analysis
 180
 181
 182   int tNormMult = 0;
 183   for (int i=0;i<nofTracks;i++)
 184     {
 185
 186       //take only primaries
 187       if(!fStack->IsPhysicalPrimary(i)) {continue;}
 188
 189       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 190
 191           //getting next track
 192       TParticle *kinetrack= fStack->Particle(i);
 193
 194       //setting multiplicity
 195         realnofTracks++;//real number of tracks (only primary particles)
 196
 197       //setting normalized multiplicity
 198       if (kinetrack->Eta() < 0.9)
 199         if(kinetrack->GetPDG()->Charge()/3!=0)
 200           tNormMult++;
 201
 202
 203           //charge
 204       trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
 205
 206
 207       //in aliroot we have AliPID
 208       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 209       //we use only 5 first
 210       double kinepid[5];
 211       for(int pid_iter=0;pid_iter<5;pid_iter++)
 212           kinepid[pid_iter]=0;
 213
 214       int pdgcode = kinetrack->GetPdgCode();
 215       //proton
 216       if(pdgcode==2212 || pdgcode==-2212)
 217         kinepid[4]=1000;
 218       //kaon
 219       if(pdgcode==321 || pdgcode==-321 )
 220         kinepid[3]=1000;
 221       //pion
 222       if( pdgcode==211 || pdgcode==-211)
 223         kinepid[2]=1000;
 224       //electron
 225       if(pdgcode==11 || pdgcode==-11)
 226         kinepid[0]=1000;
 227       //muon
 228       if(pdgcode==13 || pdgcode==-13)
 229         kinepid[1]=1000;
 230
 231       trackCopy->SetPidProbElectron(kinepid[0]);
 232       trackCopy->SetPidProbMuon(kinepid[1]);
 233       trackCopy->SetPidProbPion(kinepid[2]);
 234       trackCopy->SetPidProbKaon(kinepid[3]);
 235       trackCopy->SetPidProbProton(kinepid[4]);
 236
 237
 238         //Momentum
 239       double pxyz[3];
 240       double rxyz[3];
 241
 242         pxyz[0]=kinetrack->Px();
 243         pxyz[1]=kinetrack->Py();
 244         pxyz[2]=kinetrack->Pz();
 245
 246         rxyz[0]=kinetrack->Vx();
 247         rxyz[1]=kinetrack->Vy();
 248         rxyz[2]=kinetrack->Vz();
 249
 250         if (fRotateToEventPlane) {
 251           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 252           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 253
 254           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 255           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 256         }
 257
 258         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 259         if (v.Mag() < 0.0001) {
 260           //      cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 261           delete trackCopy;
 262           continue;
 263         }
 264
 265         trackCopy->SetP(v);//setting momentum
 266         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 267         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 268         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 269
 270         //label
 271     trackCopy->SetLabel(i);
 272
 273
 274         hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 275
 276
 277     }
 278
 279   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 280   hbtEvent->SetNormalizedMult(tNormMult);
 281   fCurEvent++;
 282
 283   return hbtEvent;
 284 }
 285
 286 //___________________
 287 void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
 288 {
 289   // The chain loads the stack for us
 290   // You must provide the address where it can be found
 291   fStack = aStack;
 292 }
 293 //___________________
 294 void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 295 {
 296   // The chain loads the generator event header for us
 297   // You must provide the address where it can be found
 298   fGenHeader = aGenHeader;
 299 }
 300
 301 //__________________
 302 void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
 303 {
 304   fRotateToEventPlane=dorotate;
 305 }
 306
 307 Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double *impact, double *covar)
 308 {
 309   // Calculates the number of sigma to the vertex.
 310
 311   Float_t b[2];
 312   Float_t bRes[2];
 313   Float_t bCov[3];
 314
 315   b[0] = impact[0];
 316   b[1] = impact[1];
 317   bCov[0] = covar[0];
 318   bCov[1] = covar[1];
 319   bCov[2] = covar[2];
 320
 321   bRes[0] = TMath::Sqrt(bCov[0]);
 322   bRes[1] = TMath::Sqrt(bCov[2]);
 323
 324   // -----------------------------------
 325   // How to get to a n-sigma cut?
 326   //
 327   // The accumulated statistics from 0 to d is
 328   //
 329   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 330   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 331   //
 332   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 333   // Can this be expressed in a different way?
 334
 335   if (bRes[0] == 0 || bRes[1] ==0)
 336     return -1;
 337
 338   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 339
 340   // stupid rounding problem screws up everything:
 341   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 342   if (TMath::Exp(-d * d / 2) < 1e-10)
 343     return 1000;
 344
 345   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 346   return d;
 347 }