[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / 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),
  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),
  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
//

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