[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / extra / AliRsnAnalysisPhi900GeV.cxx

//
// Implementation file for implementation of data analysis aft 900 GeV
//
// Author: A. Pulvirenti
//

#include "Riostream.h"

#include "TH1.h"
#include "TTree.h"
#include "TParticle.h"
#include "TRandom.h"
#include "TLorentzVector.h"

#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliStack.h"
#include "AliMCEvent.h"

#include "AliRsnAnalysisPhi900GeV.h"


AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const char *name) :
   AliAnalysisTaskSE(name),
   fUseMC(kFALSE),
   fPDG(0),
   fIM(0.0),
   fPt(0.0),
   fY(0.0),
   fEta(0.0),
   fDCAr(1E6),
   fDCAz(1E6),
   fChi2(1E6),
   fNTPC(0),
   fMinTPC(-1E6),
   fMaxTPC(1E6),
   fOutTree(),
   fOutList(0x0),
   fHEvents(0x0),
   fTOFESD(kFALSE),
   fTOFSigma(210.0),
   fTOFmaker(0x0),
   fTOFSettings(AliRsnTOFT0maker::kNone)
{
//
// Constructor
//

   DefineOutput(1, TTree::Class());
   DefineOutput(2, TTree::Class());
   DefineOutput(3, TList::Class());
}


AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const AliRsnAnalysisPhi900GeV& copy) :
   AliAnalysisTaskSE(copy),
   fUseMC(copy.fUseMC),
   fPDG(0),
   fIM(0.0),
   fPt(0.0),
   fY(0.0),
   fEta(0.0),
   fDCAr(copy.fDCAr),
   fDCAz(copy.fDCAz),
   fChi2(copy.fChi2),
   fNTPC(copy.fNTPC),
   fMinTPC(copy.fMinTPC),
   fMaxTPC(copy.fMaxTPC),
   fOutTree(),
   fOutList(0x0),
   fHEvents(0x0),
   fTOFESD(copy.fTOFESD),
   fTOFSigma(copy.fTOFSigma),
   fTOFmaker(0x0),
   fTOFSettings(copy.fTOFSettings)
{
//
// Copy constructor
//
}


AliRsnAnalysisPhi900GeV& AliRsnAnalysisPhi900GeV::operator=(const AliRsnAnalysisPhi900GeV& copy)
{
//
// Assignment operator
//
  if (this == &copy)
    return *this;
  fUseMC = copy.fUseMC;

   fDCAr = copy.fDCAr;
   fDCAz = copy.fDCAz;
   fChi2 = copy.fChi2;
   fNTPC = copy.fNTPC;

   fMinTPC = copy.fMinTPC;
   fMaxTPC = copy.fMaxTPC;

   fTOFESD      = copy.fTOFESD;
   fTOFSigma    = copy.fTOFSigma;
   fTOFSettings = copy.fTOFSettings;

   return (*this);
}


AliRsnAnalysisPhi900GeV::~AliRsnAnalysisPhi900GeV()
{
//
// Destructor
//

   if (fOutTree[0]) delete fOutTree[0];
   if (fOutTree[1]) delete fOutTree[1];
}


void AliRsnAnalysisPhi900GeV::UserCreateOutputObjects()
{
//
// Create the output data container
//

   // setup TOF maker
   fTOFmaker = new AliRsnTOFT0maker;
   fTOFmaker->SetTimeResolution(fTOFSigma * 1E-12);
   fTOFmaker->SetESDdata(fTOFESD);
   fTOFmaker->fSettings = fTOFSettings;
   AliInfo(Form("TOF sigma    = %f", fTOFSigma));
   AliInfo(Form("TOF ESD      = %s", (fTOFESD ? "YES" : "NO")));
   AliInfo(Form("TOF settings = %s", fTOFmaker->Settings().Data()));

   // load dead channel map
   fTOFmaker->LoadChannelMap((char *)"tofmap.root");
   fTOFmaker->SetMaskOffChannel();

   // initialize random
   gRandom->SetSeed(0);

   // create output trees
   OpenFile(1);
   fOutTree[0] = new TTree("rsnTree", "Pairs");

   fOutTree[0]->Branch("pdg", &fPDG, "pdg/S");
   fOutTree[0]->Branch("im" , &fIM , "im/F");
   fOutTree[0]->Branch("y"  , &fY  , "y/F");
   fOutTree[0]->Branch("pt" , &fPt , "pt/F");
   fOutTree[0]->Branch("eta", &fEta, "eta/F");

   OpenFile(2);
   fOutTree[1] = new TTree("rsnTrue", "True pairs");

   fOutTree[1]->Branch("im" , &fIM , "im/F");
   fOutTree[1]->Branch("y"  , &fY  , "y/F");
   fOutTree[1]->Branch("pt" , &fPt , "pt/F");
   fOutTree[1]->Branch("eta", &fEta, "eta/F");

   OpenFile(3);
   fOutList = new TList;
   fHEvents = new TH1I("hEvents", "Event details", 4, 0, 4);
   fOutList->Add(fHEvents);
}


void AliRsnAnalysisPhi900GeV::UserExec(Option_t *)
{
//
// Main execution function.
// Fills the fHEvents data member with the following legenda:
// 0 -- event OK, prim vertex with tracks
// 1 -- event OK, prim vertex with SPD
// 2 -- event OK but vz large
// 3 -- event bad
//

   // retrieve ESD event and related stack (if available)
   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
   AliStack    *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
   if (!esd) {
      AliError("No ESD");
      return;
   }

   // get the best primary vertex:
   // first try the one with tracks
   Int_t type = 0;
   const AliESDVertex *v = esd->GetPrimaryVertexTracks();
   if (v->GetNContributors() < 1) {
      // if not good, try SPD vertex
      type = 1;
      v = esd->GetPrimaryVertexSPD();

      // if this is not good skip this event
      if (v->GetNContributors() < 1) {
         fHEvents->Fill(3);
         PostData(3, fOutList);
         return;
      }
   }

   // if the Z position is larger than 10, skip this event
   if (TMath::Abs(v->GetZv()) > 10.0) {
      fHEvents->Fill(2);
      PostData(3, fOutList);
      return;
   }

   // use the type to fill the histogram
   fHEvents->Fill(type);

   // smear TOF times in case of MC
   if (stack) RemakeTOFtimeMC(esd);

   // get time zero for TOF
   Double_t *tof = fTOFmaker->RemakePID(esd);

   ProcessESD(esd, v, tof[0], stack);
   ProcessMC(stack);

   PostData(3, fOutList);
}


void AliRsnAnalysisPhi900GeV::Terminate(Option_t *)
{
//
// Terminate
//
}


void AliRsnAnalysisPhi900GeV::ProcessESD
(AliESDEvent *esd, const AliESDVertex *v, Double_t time0, AliStack *stack)
{
//
// This function works with the ESD object
//

   // prepare to look on all tracks to select the ones
   // which pass all the cuts
   Int_t   ntracks = esd->GetNumberOfTracks();
   TArrayI pos(ntracks);
   TArrayI neg(ntracks);

   // define fixed functions for TOF compatibility range
   Double_t a1 = 0.01, a2 = -0.03;
   Double_t b1 = 0.25, b2 =  0.25;
   Double_t c1 = 0.05, c2 = -0.03;
   Double_t ymax, ymin;

   // loop on all tracks
   Int_t    i, charge, nSPD, npos = 0, nneg = 0;
   Float_t  chi2, b[2], bCov[3];
   Double_t tpc, bb, mom, tofTime, tofRef, tofRel, times[10];
   Bool_t   okTOF;
   for (i = 0; i < ntracks; i++) {
      AliESDtrack *track = esd->GetTrack(i);
      if (!track) continue;

      // skip if it has not the required flags
      if (!track->IsOn(AliESDtrack::kTPCin)) continue;
      if (!track->IsOn(AliESDtrack::kTPCrefit)) continue;
      if (!track->IsOn(AliESDtrack::kITSrefit)) continue;

      // skip if it has not the TPC inner wall projection
      if (!track->GetInnerParam()) continue;

      // skip kink daughters
      if ((Int_t)track->GetKinkIndex(0) > 0) continue;

      // check clusters in TPC
      if (track->GetTPCclusters(0) < fNTPC) continue;

      // check chi2
      chi2  = (Float_t)track->GetTPCchi2();
      chi2 /= (Float_t)track->GetTPCclusters(0);
      if (chi2 > fChi2) continue;

      // check that has at least 1 SPD cluster
      nSPD = 0;
      if (track->HasPointOnITSLayer(0)) nSPD++;
      if (track->HasPointOnITSLayer(1)) nSPD++;
      if (nSPD < 1) continue;

      // check primary by reverting to vertex
      // and checking DCA
      if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue;
      track->GetImpactParameters(b, bCov);
      if (b[0] > fDCAr) continue;
      if (b[1] > fDCAz) continue;

      // check TPC dE/dx
      AliExternalTrackParam trackIn(*track->GetInnerParam());
      mom = trackIn.P();
      tpc = (Double_t)track->GetTPCsignal();
      bb  = AlephBB(mom);
      tpc = (tpc - bb) / bb;
      if (tpc < fMinTPC || tpc > fMaxTPC) continue;

      // if possible, check TOF
      okTOF = kTRUE;
      if (track->IsOn(AliESDtrack::kTOFpid)) {
         mom = track->P();
         if (mom <= 0.26)
            okTOF = kTRUE;
         else {
            track->GetIntegratedTimes(times);
            tofTime = (Double_t)track->GetTOFsignal() - time0;
            tofRef  = times[AliPID::kKaon];
            tofRel  = (tofTime - tofRef) / tofRef;
            ymax    = a1 / (mom - b1) + c1;
            ymin    = a2 / (mom - b2) + c2;
            okTOF   = (tofRel >= ymin && tofRel <= ymax);
         }
      }
      if (!okTOF) continue;

      // if we arrive here, all cuts were passed
      // and we add the track to one array depending on charge
      charge = (Int_t)track->Charge();
      if (charge > 0)
         pos[npos++] = i;
      else if (charge < 0)
         neg[nneg++] = i;
   }

   // resize arrays accordingly
   pos.Set(npos);
   neg.Set(nneg);

   // loop to compute invariant mass
   Int_t           ip, in, lp, ln;
   AliPID          pid;
   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
   Double_t        phimass = 1.019455;
   TParticle      *partp = 0x0, *partn = 0x0;
   AliESDtrack    *tp = 0x0, *tn = 0x0;
   TLorentzVector  vp, vn, vsum, vref;
   for (ip = 0; ip < npos; ip++) {
      tp = esd->GetTrack(pos[ip]);
      lp = TMath::Abs(tp->GetLabel());
      if (stack) partp = stack->Particle(lp);

      for (in = 0; in < nneg; in++) {
         if (pos[ip] == neg[in]) continue;
         tn = esd->GetTrack(neg[in]);
         ln = TMath::Abs(tn->GetLabel());
         if (stack) partn = stack->Particle(ln);

         fPDG = 0;
         if (partp && partn) {
            if (partp->GetFirstMother() == partn->GetFirstMother()) {
               if (partp->GetFirstMother() > 0) {
                  TParticle *mum = stack->Particle(partp->GetFirstMother());
                  fPDG = mum->GetPdgCode();
               }
            }
         }
         fPDG = TMath::Abs(fPDG);

         vp.SetXYZM(tp->Px(), tp->Py(), tp->Pz(), kmass);
         vn.SetXYZM(tn->Px(), tn->Py(), tn->Pz(), kmass);
         vsum = vp + vn;
         vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);

         fIM  = (Float_t)vsum.M();
         fPt  = (Float_t)vsum.Perp();
         fEta = (Float_t)vsum.Eta();
         fY   = (Float_t)vref.Rapidity();

         fOutTree[0]->Fill();
      }
   }

   PostData(1, fOutTree[0]);
}


void AliRsnAnalysisPhi900GeV::ProcessMC(AliStack *stack)
{
//
// Function to process stack only
//

   if (!stack) return;
   Int_t nPart = stack->GetNtrack();

   // loop to compute invariant mass
   Int_t           ip, in;
   AliPID          pid;
   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
   Double_t        phimass = 1.019455;
   TParticle      *partp = 0x0, *partn = 0x0;
   TLorentzVector  vp, vn, vsum, vref;

   for (ip = 0; ip < nPart; ip++) {
      partp = stack->Particle(ip);
      if (partp->GetPdgCode() != 321) continue;

      for (in = 0; in < nPart; in++) {
         partn = stack->Particle(in);
         if (partn->GetPdgCode() != -321) continue;

         fPDG = 0;
         if (partp->GetFirstMother() == partn->GetFirstMother()) {
            if (partp->GetFirstMother() > 0) {
               TParticle *mum = stack->Particle(partp->GetFirstMother());
               fPDG = mum->GetPdgCode();
            }
         }
         fPDG = TMath::Abs(fPDG);
         if (fPDG != 333) continue;

         vp.SetXYZM(partp->Px(), partp->Py(), partp->Pz(), kmass);
         vn.SetXYZM(partn->Px(), partn->Py(), partn->Pz(), kmass);
         vsum = vp + vn;
         vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);

         fIM  = (Float_t)vsum.M();
         fPt  = (Float_t)vsum.Perp();
         fEta = (Float_t)vsum.Eta();
         fY   = (Float_t)vref.Rapidity();

         fOutTree[1]->Fill();
      }
   }

   PostData(2, fOutTree[1]);
}


void AliRsnAnalysisPhi900GeV::SetTPCparams(Bool_t isMC)
{
//
// Set TPC bethe-bloch parameters
//

   if (!isMC) {
      fTPCpar[0] = 1.41543;
      fTPCpar[1] = 2.63394E1;
      fTPCpar[2] = 5.0411E-11;
      fTPCpar[3] = 2.12543;
      fTPCpar[4] = 4.88663;
   } else {
      fTPCpar[0] = 2.15898;
      fTPCpar[1] = 1.75295E1;
      fTPCpar[2] = 3.40030E-9;
      fTPCpar[3] = 1.96178;
      fTPCpar[4] = 3.91720;
   }
}


Double_t AliRsnAnalysisPhi900GeV::AlephBB(Double_t p, Double_t mass) const
{
//
// Compute expected Bethe-Bloch for that momentum and mass
//

   if (mass < 1E-6) return 0.0;

   Double_t aa, bb, out, beta, bg;

   bg   = p / mass;
   beta = bg / TMath::Sqrt(1.0 + bg * bg);
   aa   = TMath::Power(beta, fTPCpar[3]);
   bb   = TMath::Power(1. / bg, fTPCpar[4]);
   bb   = TMath::Log(fTPCpar[2] + bb);
   out  = (fTPCpar[1] - aa - bb) * fTPCpar[0] / aa;

   return out;
}


Double_t AliRsnAnalysisPhi900GeV::RemakeTOFtimeMC(AliESDEvent *& esd)
{
//
// Smear initial time for TOF in order to reproduce data resolution
//

   Double_t t0 = gRandom->Gaus(0, 135.); //spread in ps
   Int_t ntracks = esd->GetNumberOfTracks();
   Double_t sigmaStandard = 80.; // 80 ps from TOF

   while (ntracks--) {
      AliESDtrack *t = esd->GetTrack(ntracks);
      if ((t->GetStatus()&AliESDtrack::kTOFout) == 0 || (t->GetStatus()&AliESDtrack::kTIME) == 0) continue;
      Double_t time = t->GetTOFsignal();
      if (fTOFSigma > sigmaStandard) {
         Double_t sigmaAdded = TMath::Sqrt(fTOFSigma * fTOFSigma - sigmaStandard * sigmaStandard);
         Double_t timerandomtrack = gRandom->Gaus(0, sigmaAdded); //spread in ps
         time += timerandomtrack;
      }
      time += t0;
      t->SetTOFsignal(time);
   }
   return t0;
}
Commit	Line	Data
ca8d2905	1	//
	2	// Implementation file for implementation of data analysis aft 900 GeV
	3	//
	4	// Author: A. Pulvirenti
	5	//
	6
d7417f01	7	#include "Riostream.h"
	8
	9	#include "TH1.h"
b9bbd271	10	#include "TTree.h"
	11	#include "TParticle.h"
	12	#include "TRandom.h"
	13	#include "TLorentzVector.h"
	14
	15	#include "AliLog.h"
	16	#include "AliESDEvent.h"
	17	#include "AliESDVertex.h"
	18	#include "AliESDtrack.h"
	19	#include "AliStack.h"
	20	#include "AliMCEvent.h"
	21
	22	#include "AliRsnAnalysisPhi900GeV.h"
	23
	24
	25	AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const char *name) :
7356f978	26	AliAnalysisTaskSE(name),
	27	fUseMC(kFALSE),
	28	fPDG(0),
	29	fIM(0.0),
	30	fPt(0.0),
	31	fY(0.0),
	32	fEta(0.0),
	33	fDCAr(1E6),
	34	fDCAz(1E6),
	35	fChi2(1E6),
	36	fNTPC(0),
	37	fMinTPC(-1E6),
	38	fMaxTPC(1E6),
f5545744	39	fOutTree(),
7356f978	40	fOutList(0x0),
	41	fHEvents(0x0),
	42	fTOFESD(kFALSE),
	43	fTOFSigma(210.0),
	44	fTOFmaker(0x0),
	45	fTOFSettings(AliRsnTOFT0maker::kNone)
b9bbd271	46	{
	47	//
	48	// Constructor
	49	//
	50
7356f978	51	DefineOutput(1, TTree::Class());
	52	DefineOutput(2, TTree::Class());
	53	DefineOutput(3, TList::Class());
b9bbd271	54	}
	55
	56
	57	AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const AliRsnAnalysisPhi900GeV& copy) :
7356f978	58	AliAnalysisTaskSE(copy),
	59	fUseMC(copy.fUseMC),
	60	fPDG(0),
	61	fIM(0.0),
	62	fPt(0.0),
	63	fY(0.0),
	64	fEta(0.0),
	65	fDCAr(copy.fDCAr),
	66	fDCAz(copy.fDCAz),
	67	fChi2(copy.fChi2),
	68	fNTPC(copy.fNTPC),
	69	fMinTPC(copy.fMinTPC),
	70	fMaxTPC(copy.fMaxTPC),
f5545744	71	fOutTree(),
7356f978	72	fOutList(0x0),
	73	fHEvents(0x0),
	74	fTOFESD(copy.fTOFESD),
	75	fTOFSigma(copy.fTOFSigma),
	76	fTOFmaker(0x0),
	77	fTOFSettings(copy.fTOFSettings)
b9bbd271	78	{
	79	//
	80	// Copy constructor
	81	//
	82	}
	83
	84
	85	AliRsnAnalysisPhi900GeV& AliRsnAnalysisPhi900GeV::operator=(const AliRsnAnalysisPhi900GeV& copy)
	86	{
	87	//
	88	// Assignment operator
	89	//
aaf7af8d	90	if (this == &copy)
	91	return *this;
	92	fUseMC = copy.fUseMC;
b9bbd271	93
7356f978	94	fDCAr = copy.fDCAr;
	95	fDCAz = copy.fDCAz;
	96	fChi2 = copy.fChi2;
	97	fNTPC = copy.fNTPC;
b9bbd271	98
7356f978	99	fMinTPC = copy.fMinTPC;
7356f978	100	fMaxTPC = copy.fMaxTPC;
b9bbd271	101
7356f978	102	fTOFESD = copy.fTOFESD;
	103	fTOFSigma = copy.fTOFSigma;
	104	fTOFSettings = copy.fTOFSettings;
b9bbd271	105
7356f978	106	return (*this);
b9bbd271	107	}
	108
	109
	110	AliRsnAnalysisPhi900GeV::~AliRsnAnalysisPhi900GeV()
	111	{
	112	//
	113	// Destructor
	114	//
	115
7356f978	116	if (fOutTree[0]) delete fOutTree[0];
7356f978	117	if (fOutTree[1]) delete fOutTree[1];
b9bbd271	118	}
	119
	120
	121	void AliRsnAnalysisPhi900GeV::UserCreateOutputObjects()
	122	{
	123	//
	124	// Create the output data container
	125	//
	126
7356f978	127	// setup TOF maker
	128	fTOFmaker = new AliRsnTOFT0maker;
	129	fTOFmaker->SetTimeResolution(fTOFSigma * 1E-12);
	130	fTOFmaker->SetESDdata(fTOFESD);
	131	fTOFmaker->fSettings = fTOFSettings;
	132	AliInfo(Form("TOF sigma = %f", fTOFSigma));
	133	AliInfo(Form("TOF ESD = %s", (fTOFESD ? "YES" : "NO")));
	134	AliInfo(Form("TOF settings = %s", fTOFmaker->Settings().Data()));
	135
	136	// load dead channel map
	137	fTOFmaker->LoadChannelMap((char *)"tofmap.root");
	138	fTOFmaker->SetMaskOffChannel();
	139
	140	// initialize random
	141	gRandom->SetSeed(0);
	142
	143	// create output trees
	144	OpenFile(1);
	145	fOutTree[0] = new TTree("rsnTree", "Pairs");
	146
	147	fOutTree[0]->Branch("pdg", &fPDG, "pdg/S");
	148	fOutTree[0]->Branch("im" , &fIM , "im/F");
	149	fOutTree[0]->Branch("y" , &fY , "y/F");
	150	fOutTree[0]->Branch("pt" , &fPt , "pt/F");
	151	fOutTree[0]->Branch("eta", &fEta, "eta/F");
	152
	153	OpenFile(2);
	154	fOutTree[1] = new TTree("rsnTrue", "True pairs");
	155
	156	fOutTree[1]->Branch("im" , &fIM , "im/F");
	157	fOutTree[1]->Branch("y" , &fY , "y/F");
	158	fOutTree[1]->Branch("pt" , &fPt , "pt/F");
	159	fOutTree[1]->Branch("eta", &fEta, "eta/F");
	160
	161	OpenFile(3);
	162	fOutList = new TList;
	163	fHEvents = new TH1I("hEvents", "Event details", 4, 0, 4);
	164	fOutList->Add(fHEvents);
b9bbd271	165	}
	166
	167
	168	void AliRsnAnalysisPhi900GeV::UserExec(Option_t *)
	169	{
	170	//
d7417f01	171	// Main execution function.
	172	// Fills the fHEvents data member with the following legenda:
	173	// 0 -- event OK, prim vertex with tracks
	174	// 1 -- event OK, prim vertex with SPD
	175	// 2 -- event OK but vz large
	176	// 3 -- event bad
b9bbd271	177	//
b9bbd271	178
7356f978	179	// retrieve ESD event and related stack (if available)
	180	AliESDEvent esd = dynamic_cast<AliESDEvent>(fInputEvent);
	181	AliStack *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
245cd328	182	if (!esd) {
	183	AliError("No ESD");
	184	return;
	185	}
7356f978	186
	187	// get the best primary vertex:
	188	// first try the one with tracks
	189	Int_t type = 0;
	190	const AliESDVertex *v = esd->GetPrimaryVertexTracks();
	191	if (v->GetNContributors() < 1) {
	192	// if not good, try SPD vertex
	193	type = 1;
	194	v = esd->GetPrimaryVertexSPD();
	195
	196	// if this is not good skip this event
	197	if (v->GetNContributors() < 1) {
	198	fHEvents->Fill(3);
	199	PostData(3, fOutList);
	200	return;
	201	}
	202	}
	203
	204	// if the Z position is larger than 10, skip this event
	205	if (TMath::Abs(v->GetZv()) > 10.0) {
	206	fHEvents->Fill(2);
d7417f01	207	PostData(3, fOutList);
d7417f01	208	return;
7356f978	209	}
d7417f01	210
7356f978	211	// use the type to fill the histogram
7356f978	212	fHEvents->Fill(type);
b9bbd271	213
7356f978	214	// smear TOF times in case of MC
7356f978	215	if (stack) RemakeTOFtimeMC(esd);
b9bbd271	216
7356f978	217	// get time zero for TOF
7356f978	218	Double_t *tof = fTOFmaker->RemakePID(esd);
b9bbd271	219
7356f978	220	ProcessESD(esd, v, tof[0], stack);
7356f978	221	ProcessMC(stack);
d7417f01	222
7356f978	223	PostData(3, fOutList);
b9bbd271	224	}
	225
	226
	227	void AliRsnAnalysisPhi900GeV::Terminate(Option_t *)
	228	{
	229	//
	230	// Terminate
	231	//
	232	}
	233
	234
	235	void AliRsnAnalysisPhi900GeV::ProcessESD
	236	(AliESDEvent esd, const AliESDVertex v, Double_t time0, AliStack *stack)
	237	{
	238	//
	239	// This function works with the ESD object
	240	//
	241
7356f978	242	// prepare to look on all tracks to select the ones
	243	// which pass all the cuts
	244	Int_t ntracks = esd->GetNumberOfTracks();
	245	TArrayI pos(ntracks);
	246	TArrayI neg(ntracks);
	247
	248	// define fixed functions for TOF compatibility range
	249	Double_t a1 = 0.01, a2 = -0.03;
	250	Double_t b1 = 0.25, b2 = 0.25;
	251	Double_t c1 = 0.05, c2 = -0.03;
	252	Double_t ymax, ymin;
	253
	254	// loop on all tracks
	255	Int_t i, charge, nSPD, npos = 0, nneg = 0;
	256	Float_t chi2, b[2], bCov[3];
	257	Double_t tpc, bb, mom, tofTime, tofRef, tofRel, times[10];
	258	Bool_t okTOF;
	259	for (i = 0; i < ntracks; i++) {
	260	AliESDtrack *track = esd->GetTrack(i);
	261	if (!track) continue;
	262
	263	// skip if it has not the required flags
	264	if (!track->IsOn(AliESDtrack::kTPCin)) continue;
	265	if (!track->IsOn(AliESDtrack::kTPCrefit)) continue;
	266	if (!track->IsOn(AliESDtrack::kITSrefit)) continue;
	267
	268	// skip if it has not the TPC inner wall projection
	269	if (!track->GetInnerParam()) continue;
	270
	271	// skip kink daughters
	272	if ((Int_t)track->GetKinkIndex(0) > 0) continue;
	273
	274	// check clusters in TPC
	275	if (track->GetTPCclusters(0) < fNTPC) continue;
	276
	277	// check chi2
	278	chi2 = (Float_t)track->GetTPCchi2();
	279	chi2 /= (Float_t)track->GetTPCclusters(0);
	280	if (chi2 > fChi2) continue;
	281
	282	// check that has at least 1 SPD cluster
	283	nSPD = 0;
	284	if (track->HasPointOnITSLayer(0)) nSPD++;
	285	if (track->HasPointOnITSLayer(1)) nSPD++;
	286	if (nSPD < 1) continue;
	287
	288	// check primary by reverting to vertex
	289	// and checking DCA
	290	if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue;
	291	track->GetImpactParameters(b, bCov);
	292	if (b[0] > fDCAr) continue;
	293	if (b[1] > fDCAz) continue;
	294
	295	// check TPC dE/dx
	296	AliExternalTrackParam trackIn(*track->GetInnerParam());
	297	mom = trackIn.P();
	298	tpc = (Double_t)track->GetTPCsignal();
	299	bb = AlephBB(mom);
	300	tpc = (tpc - bb) / bb;
	301	if (tpc < fMinTPC \|\| tpc > fMaxTPC) continue;
	302
	303	// if possible, check TOF
	304	okTOF = kTRUE;
	305	if (track->IsOn(AliESDtrack::kTOFpid)) {
306	mom = track->P();
307	if (mom <= 0.26)
308	okTOF = kTRUE;
309	else {
310	track->GetIntegratedTimes(times);
311	tofTime = (Double_t)track->GetTOFsignal() - time0;
312	tofRef = times[AliPID::kKaon];
313	tofRel = (tofTime - tofRef) / tofRef;
314	ymax = a1 / (mom - b1) + c1;
315	ymin = a2 / (mom - b2) + c2;
316	okTOF = (tofRel >= ymin && tofRel <= ymax);
317	}
b9bbd271	318	}
7356f978	319	if (!okTOF) continue;
	320
	321	// if we arrive here, all cuts were passed
	322	// and we add the track to one array depending on charge
	323	charge = (Int_t)track->Charge();
	324	if (charge > 0)
	325	pos[npos++] = i;
	326	else if (charge < 0)
	327	neg[nneg++] = i;
	328	}
b9bbd271	329
7356f978	330	// resize arrays accordingly
	331	pos.Set(npos);
	332	neg.Set(nneg);
	333
	334	// loop to compute invariant mass
	335	Int_t ip, in, lp, ln;
	336	AliPID pid;
	337	Double_t kmass = pid.ParticleMass(AliPID::kKaon);
	338	Double_t phimass = 1.019455;
	339	TParticle partp = 0x0, partn = 0x0;
	340	AliESDtrack tp = 0x0, tn = 0x0;
	341	TLorentzVector vp, vn, vsum, vref;
	342	for (ip = 0; ip < npos; ip++) {
	343	tp = esd->GetTrack(pos[ip]);
	344	lp = TMath::Abs(tp->GetLabel());
	345	if (stack) partp = stack->Particle(lp);
	346
	347	for (in = 0; in < nneg; in++) {
	348	if (pos[ip] == neg[in]) continue;
	349	tn = esd->GetTrack(neg[in]);
	350	ln = TMath::Abs(tn->GetLabel());
	351	if (stack) partn = stack->Particle(ln);
	352
	353	fPDG = 0;
	354	if (partp && partn) {
	355	if (partp->GetFirstMother() == partn->GetFirstMother()) {
	356	if (partp->GetFirstMother() > 0) {
	357	TParticle *mum = stack->Particle(partp->GetFirstMother());
	358	fPDG = mum->GetPdgCode();
	359	}
	360	}
	361	}
	362	fPDG = TMath::Abs(fPDG);
	363
	364	vp.SetXYZM(tp->Px(), tp->Py(), tp->Pz(), kmass);
	365	vn.SetXYZM(tn->Px(), tn->Py(), tn->Pz(), kmass);
	366	vsum = vp + vn;
	367	vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
	368
	369	fIM = (Float_t)vsum.M();
	370	fPt = (Float_t)vsum.Perp();
	371	fEta = (Float_t)vsum.Eta();
	372	fY = (Float_t)vref.Rapidity();
	373
	374	fOutTree[0]->Fill();
	375	}
	376	}
d7417f01	377
7356f978	378	PostData(1, fOutTree[0]);
b9bbd271	379	}
	380
	381
	382	void AliRsnAnalysisPhi900GeV::ProcessMC(AliStack *stack)
	383	{
	384	//
	385	// Function to process stack only
	386	//
	387
7356f978	388	if (!stack) return;
	389	Int_t nPart = stack->GetNtrack();
	390
	391	// loop to compute invariant mass
	392	Int_t ip, in;
	393	AliPID pid;
	394	Double_t kmass = pid.ParticleMass(AliPID::kKaon);
	395	Double_t phimass = 1.019455;
	396	TParticle partp = 0x0, partn = 0x0;
	397	TLorentzVector vp, vn, vsum, vref;
	398
	399	for (ip = 0; ip < nPart; ip++) {
	400	partp = stack->Particle(ip);
	401	if (partp->GetPdgCode() != 321) continue;
	402
	403	for (in = 0; in < nPart; in++) {
	404	partn = stack->Particle(in);
	405	if (partn->GetPdgCode() != -321) continue;
	406
	407	fPDG = 0;
	408	if (partp->GetFirstMother() == partn->GetFirstMother()) {
	409	if (partp->GetFirstMother() > 0) {
	410	TParticle *mum = stack->Particle(partp->GetFirstMother());
	411	fPDG = mum->GetPdgCode();
	412	}
	413	}
	414	fPDG = TMath::Abs(fPDG);
	415	if (fPDG != 333) continue;
	416
	417	vp.SetXYZM(partp->Px(), partp->Py(), partp->Pz(), kmass);
	418	vn.SetXYZM(partn->Px(), partn->Py(), partn->Pz(), kmass);
	419	vsum = vp + vn;
	420	vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
	421
	422	fIM = (Float_t)vsum.M();
	423	fPt = (Float_t)vsum.Perp();
	424	fEta = (Float_t)vsum.Eta();
	425	fY = (Float_t)vref.Rapidity();
	426
	427	fOutTree[1]->Fill();
b9bbd271	428	}
7356f978	429	}
d7417f01	430
7356f978	431	PostData(2, fOutTree[1]);
d7417f01	432	}
	433
	434
	435	void AliRsnAnalysisPhi900GeV::SetTPCparams(Bool_t isMC)
	436	{
	437	//
	438	// Set TPC bethe-bloch parameters
	439	//
	440
7356f978	441	if (!isMC) {
	442	fTPCpar[0] = 1.41543;
	443	fTPCpar[1] = 2.63394E1;
	444	fTPCpar[2] = 5.0411E-11;
	445	fTPCpar[3] = 2.12543;
	446	fTPCpar[4] = 4.88663;
	447	} else {
	448	fTPCpar[0] = 2.15898;
	449	fTPCpar[1] = 1.75295E1;
	450	fTPCpar[2] = 3.40030E-9;
	451	fTPCpar[3] = 1.96178;
	452	fTPCpar[4] = 3.91720;
	453	}
b9bbd271	454	}
	455
	456
ca8d2905	457	Double_t AliRsnAnalysisPhi900GeV::AlephBB(Double_t p, Double_t mass) const
b9bbd271	458	{
	459	//
	460	// Compute expected Bethe-Bloch for that momentum and mass
	461	//
	462
7356f978	463	if (mass < 1E-6) return 0.0;
	464
	465	Double_t aa, bb, out, beta, bg;
b9bbd271	466
7356f978	467	bg = p / mass;
	468	beta = bg / TMath::Sqrt(1.0 + bg * bg);
	469	aa = TMath::Power(beta, fTPCpar[3]);
	470	bb = TMath::Power(1. / bg, fTPCpar[4]);
	471	bb = TMath::Log(fTPCpar[2] + bb);
	472	out = (fTPCpar[1] - aa - bb) * fTPCpar[0] / aa;
b9bbd271	473
7356f978	474	return out;
b9bbd271	475	}
	476
	477
	478	Double_t AliRsnAnalysisPhi900GeV::RemakeTOFtimeMC(AliESDEvent *& esd)
	479	{
	480	//
	481	// Smear initial time for TOF in order to reproduce data resolution
	482	//
	483
7356f978	484	Double_t t0 = gRandom->Gaus(0, 135.); //spread in ps
	485	Int_t ntracks = esd->GetNumberOfTracks();
	486	Double_t sigmaStandard = 80.; // 80 ps from TOF
	487
	488	while (ntracks--) {
	489	AliESDtrack *t = esd->GetTrack(ntracks);
	490	if ((t->GetStatus()&AliESDtrack::kTOFout) == 0 \|\| (t->GetStatus()&AliESDtrack::kTIME) == 0) continue;
	491	Double_t time = t->GetTOFsignal();
	492	if (fTOFSigma > sigmaStandard) {
	493	Double_t sigmaAdded = TMath::Sqrt(fTOFSigma * fTOFSigma - sigmaStandard * sigmaStandard);
	494	Double_t timerandomtrack = gRandom->Gaus(0, sigmaAdded); //spread in ps
	495	time += timerandomtrack;
	496	}
	497	time += t0;
	498	t->SetTOFsignal(time);
b9bbd271	499	}
7356f978	500	return t0;
7356f978	501	}