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