[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetFastSimulation.cxx

// $Id$
//
// Jet model task to merge to existing branches
// only implemented for track branches
//
// Author: M. Verweij

#include "AliJetFastSimulation.h"

#include <TClonesArray.h>
#include <TFolder.h>
#include <TLorentzVector.h>
#include <TParticle.h>
#include <TParticlePDG.h>
#include <TRandom3.h>
#include <TProfile.h>
#include <TGrid.h>
#include <TFile.h>
#include <TF1.h>
#include "AliAnalysisManager.h"
#include "AliEMCALDigit.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALRecPoint.h"
#include "AliGenerator.h"
#include "AliHeader.h"
#include "AliLog.h"
#include "AliPicoTrack.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliStack.h"
#include "AliStack.h"
#include "AliVCluster.h"
#include "AliVEvent.h"

ClassImp(AliJetFastSimulation)

//________________________________________________________________________
AliJetFastSimulation::AliJetFastSimulation() : 
AliAnalysisTaskEmcal("AliJetFastSimulation",kTRUE),
  fTracksOutName(""),
  fTracksOut(0x0),
  fNTrackClasses(2),
  fRandom(0),
  fEfficiencyFixed(1.1),
  fMomResH1(0x0),
  fMomResH2(0x0),
  fMomResH3(0x0),
  fMomResH1Fit(0x0),
  fMomResH2Fit(0x0),
  fMomResH3Fit(0x0),
  fhEffH1(0x0),
  fhEffH2(0x0),
  fhEffH3(0x0),
  fUseTrPtResolutionSmearing(kFALSE),
  fUseDiceEfficiency(kFALSE),
  fDiceEfficiencyMinPt(-1.),
  fUseTrPtResolutionFromOADB(kFALSE),
  fUseTrEfficiencyFromOADB(kFALSE),
  fPathTrPtResolution(""),
  fPathTrEfficiency(""),
  fHistPtDet(0),
  fh2PtGenPtSmeared(0),
  fp1Efficiency(0),
  fp1PtResolution(0)
{
  // Default constructor.
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliJetFastSimulation::AliJetFastSimulation(const char *name) : 
  AliAnalysisTaskEmcal(name,kTRUE),
  fTracksOutName(""),
  fTracksOut(0x0),
  fNTrackClasses(2),
  fRandom(0),
  fEfficiencyFixed(1.1),
  fMomResH1(0x0),
  fMomResH2(0x0),
  fMomResH3(0x0),
  fMomResH1Fit(0x0),
  fMomResH2Fit(0x0),
  fMomResH3Fit(0x0),
  fhEffH1(0x0),
  fhEffH2(0x0),
  fhEffH3(0x0),
  fUseTrPtResolutionSmearing(kFALSE),
  fUseDiceEfficiency(kFALSE),
  fDiceEfficiencyMinPt(-1.),
  fUseTrPtResolutionFromOADB(kFALSE),
  fUseTrEfficiencyFromOADB(kFALSE),
  fPathTrPtResolution(""),
  fPathTrEfficiency(""),
  fHistPtDet(0),
  fh2PtGenPtSmeared(0),
  fp1Efficiency(0),
  fp1PtResolution(0)
{
  // Standard constructor.
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliJetFastSimulation::~AliJetFastSimulation()
{
  // Destructor

  delete fRandom;
}

//________________________________________________________________________
void AliJetFastSimulation::ExecOnce() 
{
  // Exec only once.

  AliAnalysisTaskEmcal::ExecOnce();

  if(!fRandom) fRandom = new TRandom3(0);

  if (!fTracksOutName.IsNull()) {
    fTracksOut = new TClonesArray("AliPicoTrack");
    fTracksOut->SetName(fTracksOutName);
    if (InputEvent()->FindListObject(fTracksOutName)) {
      AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fTracksOutName.Data()));
      return;
    }
    else {
      InputEvent()->AddObject(fTracksOut);
    }
  }
}
//________________________________________________________________________
void AliJetFastSimulation::LocalInit() {
  //initialize track response
  if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
  if(fUseTrEfficiencyFromOADB)   LoadTrEfficiencyRootFileFromOADB();
  FitMomentumResolution();
}

//________________________________________________________________________
void AliJetFastSimulation::UserCreateOutputObjects() 
{
  AliAnalysisTaskEmcal::UserCreateOutputObjects();

  const Int_t nBinPt = 100;
  Double_t binLimitsPt[nBinPt+1];
  for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
    if(iPt == 0){
      binLimitsPt[iPt] = 0.0;
    } else {// 1.0
      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
    }
  }

  fHistPtDet = new TH1F("fHistpt","fHistPtDet;#it{p}_{T};N",nBinPt,binLimitsPt);
  fOutput->Add(fHistPtDet);

  fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
  fOutput->Add(fh2PtGenPtSmeared);

  fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
  fOutput->Add(fp1Efficiency);

  fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
  fOutput->Add(fp1PtResolution);

  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
}


//________________________________________________________________________
Bool_t AliJetFastSimulation::Run() 
{
  //Check if information is provided detector level effects
  if(!fMomResH1 && fNTrackClasses>0 ) 
    fUseTrPtResolutionSmearing = 0;
  if(!fMomResH2 && fNTrackClasses>1 ) 
    fUseTrPtResolutionSmearing = 0;
  if(!fMomResH3 && fNTrackClasses>2 ) 
    fUseTrPtResolutionSmearing = 0;

  if(fEfficiencyFixed < 1.)
    fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user but not implemented
  else {
    if(!fhEffH1 && fNTrackClasses>0 ) 
      fUseDiceEfficiency = 0;
    if(!fhEffH2 && fNTrackClasses>1 ) 
      fUseDiceEfficiency = 0;
    if(!fhEffH3 && fNTrackClasses>2 ) 
      fUseDiceEfficiency = 0;
  }

  fTracksOut->Delete();
  SimulateTracks();
  return kTRUE;
}

//________________________________________________________________________
void AliJetFastSimulation::SimulateTracks()
{
  //Apply toy detector simulation to tracks
  Int_t it = 0;
  const Int_t nTracks = fTracks->GetEntriesFast();
   for (Int_t i = 0; i < nTracks; ++i) {
    AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
    if (!picotrack)
      continue;

    Bool_t accept = kTRUE;
    Double_t eff[3] = {0};
    Double_t rnd = fRandom->Uniform(1.);  
    if(fUseDiceEfficiency) accept = DiceEfficiency(picotrack,eff,rnd);
    if(!accept) continue;

    AliPicoTrack *track = NULL;
    if(fUseTrPtResolutionSmearing) {
      track = SmearPt(picotrack,eff,rnd);
      (*fTracksOut)[it] = track;
    } else
      track = new ((*fTracksOut)[it]) AliPicoTrack(*picotrack);

    track->SetBit(TObject::kBitMask,1);
    fHistPtDet->Fill(track->Pt());
    it++;
   }
}

//________________________________________________________________________
Bool_t AliJetFastSimulation::DiceEfficiency(AliPicoTrack *vp, Double_t eff[3], Double_t rnd)
{
  // Dice to decide if particle is kept or not - toy  model for efficiency
  //
  Double_t sumEff = 0.;
  Double_t pT = 0.;
 
  if(fEfficiencyFixed<1.)
    sumEff = fEfficiencyFixed;
  else {
    pT = vp->Pt();
    Double_t pTtmp = pT;
    if(pT>10.) pTtmp = 10.;
    if(fhEffH1) eff[0] = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
    if(fhEffH2) eff[1] = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
    if(fhEffH3) eff[2] = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));

    sumEff = eff[0]+eff[1]+eff[2];
  }
  fp1Efficiency->Fill(vp->Pt(),sumEff);
  if(rnd>sumEff && pT > fDiceEfficiencyMinPt) return kFALSE;
  return kTRUE;
}

//________________________________________________________________________
AliPicoTrack* AliJetFastSimulation::SmearPt(AliPicoTrack *vp, Double_t eff[3], Double_t rnd)
{
  //Smear momentum of generated particle
  Double_t smear = 1.;
  Double_t  pT = vp->Pt();
  //Select hybrid track category

  //Sort efficiencies from large to small
  Int_t cat[3] = {0};
  TMath::Sort(3,eff,cat);
  if(rnd<=eff[cat[2]])
    smear = GetMomentumSmearing(cat[2],pT);
  else if(rnd<=(eff[cat[2]]+eff[cat[1]]))
    smear = GetMomentumSmearing(cat[1],pT);
  else
    smear = GetMomentumSmearing(cat[0],pT);
  
  fp1PtResolution->Fill(vp->Pt(),smear);

  Double_t sigma = vp->Pt()*smear;
  Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
  fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);

  AliPicoTrack *picotrack = new AliPicoTrack(pTrec,
					     vp->Eta(),
					     vp->Phi(),
					     vp->Charge(),
					     vp->GetLabel(),
					     AliPicoTrack::GetTrackType(vp),
					     vp->GetTrackEtaOnEMCal(),
					     vp->GetTrackPhiOnEMCal(),
					     vp->GetTrackPtOnEMCal(),
					     vp->IsEMCAL(),
					     0.13957); //assume pion mass
    return picotrack;
}

//________________________________________________________________________
Double_t AliJetFastSimulation::GetMomentumSmearing(Int_t cat, Double_t pt) {

  //
  // Get smearing on generated momentum
  //

  TProfile *fMomRes = 0x0;
  if(cat==1 && fMomResH1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
  if(cat==2 && fMomResH2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
  if(cat==3 && fMomResH3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");

  if(!fMomRes)
    return 0.;

  Double_t smear = 0.;
  if(pt>20.) {
    if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
    if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
    if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
  }
  else {
    Int_t bin = fMomRes->FindBin(pt);
    smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
  }
  
  if(fMomRes) delete fMomRes;

  return smear;
}

//________________________________________________________________________
void AliJetFastSimulation::LoadTrPtResolutionRootFileFromOADB() {

  if (!gGrid && fPathTrPtResolution.Contains("alien://")) {
     AliInfo("Trying to connect to AliEn ...");
     TGrid::Connect("alien://");
   }

  TFile *f = TFile::Open(fPathTrPtResolution.Data());
  if(!f)return;
  TProfile *fProfPtPtSigma1PtGlobSt     = NULL;
  TProfile *fProfPtPtSigma1PtGlobCnoSPD = NULL;
  TProfile *fProfPtPtSigma1PtGlobCnoITS = NULL;
  if(fNTrackClasses>0) fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
  if(fNTrackClasses>1) fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
  if(fNTrackClasses>2) fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");

  SetSmearResolution(kTRUE);
  SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoSPD,fProfPtPtSigma1PtGlobCnoITS);
}

//________________________________________________________________________
void AliJetFastSimulation::LoadTrEfficiencyRootFileFromOADB() {

   if (!gGrid && fPathTrPtResolution.Contains("alien://")) {
     AliInfo("Trying to connect to AliEn ...");
     TGrid::Connect("alien://");
   }

  TFile *f = TFile::Open(fPathTrEfficiency.Data());
  if(!f)return;
  TH1D *hEffPosGlobSt     = NULL;
  TH1D *hEffPosGlobCnoSPD = NULL;
  TH1D *hEffPosGlobCnoITS = NULL;
  if(fNTrackClasses>0) hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
  if(fNTrackClasses>1) hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
  if(fNTrackClasses>2) hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");

  SetDiceEfficiency(kTRUE);
  SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoSPD,hEffPosGlobCnoITS);
}

//________________________________________________________________________
void AliJetFastSimulation::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
  //
  // set mom res profiles
  //
  if(fMomResH1) delete fMomResH1;
  if(fMomResH2) delete fMomResH2;
  if(fMomResH3) delete fMomResH3;
  
  if(p1) fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
  if(p2) fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
  if(p3) fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
}

//________________________________________________________________________
void AliJetFastSimulation:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
  //
  // set tracking efficiency histos
  //
  if(h1) fhEffH1 = (TH1*)h1->Clone("fhEffH1");
  if(h2) fhEffH2 = (TH1*)h2->Clone("fhEffH2");
  if(h3) fhEffH3 = (TH1*)h3->Clone("fhEffH3");
}

//________________________________________________________________________
void AliJetFastSimulation::FitMomentumResolution() {
  //
  // Fit linear function on momentum resolution at high pT
  //

  if(!fMomResH1Fit && fMomResH1) {
    fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
    fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
    fMomResH1Fit ->SetRange(5.,100.);
  }

  if(!fMomResH2Fit && fMomResH2) {
    fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
    fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
    fMomResH2Fit ->SetRange(5.,100.);
  }

  if(!fMomResH3Fit && fMomResH3) {
    fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
    fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
    fMomResH3Fit ->SetRange(5.,100.);
  }

}
Commit	Line	Data
b3821b05	1	// $Id$
	2	//
	3	// Jet model task to merge to existing branches
	4	// only implemented for track branches
	5	//
	6	// Author: M. Verweij
	7
	8	#include "AliJetFastSimulation.h"
	9
	10	#include <TClonesArray.h>
	11	#include <TFolder.h>
	12	#include <TLorentzVector.h>
	13	#include <TParticle.h>
	14	#include <TParticlePDG.h>
	15	#include <TRandom3.h>
	16	#include <TProfile.h>
	17	#include <TGrid.h>
	18	#include <TFile.h>
	19	#include <TF1.h>
	20	#include "AliAnalysisManager.h"
	21	#include "AliEMCALDigit.h"
	22	#include "AliEMCALGeometry.h"
	23	#include "AliEMCALRecPoint.h"
	24	#include "AliGenerator.h"
	25	#include "AliHeader.h"
	26	#include "AliLog.h"
	27	#include "AliPicoTrack.h"
	28	#include "AliRun.h"
	29	#include "AliRunLoader.h"
	30	#include "AliStack.h"
	31	#include "AliStack.h"
	32	#include "AliVCluster.h"
	33	#include "AliVEvent.h"
	34
	35	ClassImp(AliJetFastSimulation)
	36
	37	//________________________________________________________________________
	38	AliJetFastSimulation::AliJetFastSimulation() :
	39	AliAnalysisTaskEmcal("AliJetFastSimulation",kTRUE),
	40	fTracksOutName(""),
	41	fTracksOut(0x0),
	42	fNTrackClasses(2),
	43	fRandom(0),
689017c8	44	fEfficiencyFixed(1.1),
b3821b05	45	fMomResH1(0x0),
	46	fMomResH2(0x0),
	47	fMomResH3(0x0),
	48	fMomResH1Fit(0x0),
	49	fMomResH2Fit(0x0),
	50	fMomResH3Fit(0x0),
	51	fhEffH1(0x0),
	52	fhEffH2(0x0),
	53	fhEffH3(0x0),
	54	fUseTrPtResolutionSmearing(kFALSE),
	55	fUseDiceEfficiency(kFALSE),
	56	fDiceEfficiencyMinPt(-1.),
	57	fUseTrPtResolutionFromOADB(kFALSE),
	58	fUseTrEfficiencyFromOADB(kFALSE),
	59	fPathTrPtResolution(""),
	60	fPathTrEfficiency(""),
	61	fHistPtDet(0),
	62	fh2PtGenPtSmeared(0),
	63	fp1Efficiency(0),
	64	fp1PtResolution(0)
	65	{
	66	// Default constructor.
	67	SetMakeGeneralHistograms(kTRUE);
	68	}
	69
	70	//________________________________________________________________________
	71	AliJetFastSimulation::AliJetFastSimulation(const char *name) :
	72	AliAnalysisTaskEmcal(name,kTRUE),
	73	fTracksOutName(""),
	74	fTracksOut(0x0),
	75	fNTrackClasses(2),
	76	fRandom(0),
689017c8	77	fEfficiencyFixed(1.1),
b3821b05	78	fMomResH1(0x0),
	79	fMomResH2(0x0),
	80	fMomResH3(0x0),
	81	fMomResH1Fit(0x0),
	82	fMomResH2Fit(0x0),
	83	fMomResH3Fit(0x0),
	84	fhEffH1(0x0),
	85	fhEffH2(0x0),
	86	fhEffH3(0x0),
	87	fUseTrPtResolutionSmearing(kFALSE),
	88	fUseDiceEfficiency(kFALSE),
	89	fDiceEfficiencyMinPt(-1.),
	90	fUseTrPtResolutionFromOADB(kFALSE),
	91	fUseTrEfficiencyFromOADB(kFALSE),
	92	fPathTrPtResolution(""),
	93	fPathTrEfficiency(""),
	94	fHistPtDet(0),
	95	fh2PtGenPtSmeared(0),
	96	fp1Efficiency(0),
	97	fp1PtResolution(0)
	98	{
	99	// Standard constructor.
	100	SetMakeGeneralHistograms(kTRUE);
	101	}
	102
	103	//________________________________________________________________________
	104	AliJetFastSimulation::~AliJetFastSimulation()
	105	{
	106	// Destructor
	107
	108	delete fRandom;
	109	}
	110
	111	//________________________________________________________________________
	112	void AliJetFastSimulation::ExecOnce()
	113	{
	114	// Exec only once.
	115
	116	AliAnalysisTaskEmcal::ExecOnce();
	117
	118	if(!fRandom) fRandom = new TRandom3(0);
	119
	120	if (!fTracksOutName.IsNull()) {
	121	fTracksOut = new TClonesArray("AliPicoTrack");
	122	fTracksOut->SetName(fTracksOutName);
	123	if (InputEvent()->FindListObject(fTracksOutName)) {
	124	AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fTracksOutName.Data()));
	125	return;
	126	}
	127	else {
	128	InputEvent()->AddObject(fTracksOut);
	129	}
	130	}
	131	}
	132	//________________________________________________________________________
	133	void AliJetFastSimulation::LocalInit() {
	134	//initialize track response
	135	if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
	136	if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
	137	FitMomentumResolution();
	138	}
	139
	140	//________________________________________________________________________
	141	void AliJetFastSimulation::UserCreateOutputObjects()
142	{
143	AliAnalysisTaskEmcal::UserCreateOutputObjects();
144
145	const Int_t nBinPt = 100;
146	Double_t binLimitsPt[nBinPt+1];
147	for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
148	if(iPt == 0){
149	binLimitsPt[iPt] = 0.0;
150	} else {// 1.0
151	binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
152	}
153	}
154
155	fHistPtDet = new TH1F("fHistpt","fHistPtDet;#it{p}_{T};N",nBinPt,binLimitsPt);
156	fOutput->Add(fHistPtDet);
157
158	fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
159	fOutput->Add(fh2PtGenPtSmeared);
160
161	fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
162	fOutput->Add(fp1Efficiency);
163
164	fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
165	fOutput->Add(fp1PtResolution);
166
167	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
168	}
169
170
171	//________________________________________________________________________
172	Bool_t AliJetFastSimulation::Run()
173	{
174	//Check if information is provided detector level effects
175	if(!fMomResH1 && fNTrackClasses>0 )
176	fUseTrPtResolutionSmearing = 0;
177	if(!fMomResH2 && fNTrackClasses>1 )
178	fUseTrPtResolutionSmearing = 0;
179	if(!fMomResH3 && fNTrackClasses>2 )
180	fUseTrPtResolutionSmearing = 0;
181
689017c8	182	if(fEfficiencyFixed < 1.)
b3821b05	183	fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user but not implemented
	184	else {
	185	if(!fhEffH1 && fNTrackClasses>0 )
	186	fUseDiceEfficiency = 0;
	187	if(!fhEffH2 && fNTrackClasses>1 )
	188	fUseDiceEfficiency = 0;
	189	if(!fhEffH3 && fNTrackClasses>2 )
	190	fUseDiceEfficiency = 0;
	191	}
	192
3c923384	193	fTracksOut->Delete();
b3821b05	194	SimulateTracks();
	195	return kTRUE;
	196	}
	197
	198	//________________________________________________________________________
	199	void AliJetFastSimulation::SimulateTracks()
	200	{
	201	//Apply toy detector simulation to tracks
4ccfe90c	202	Int_t it = 0;
b3821b05	203	const Int_t nTracks = fTracks->GetEntriesFast();
	204	for (Int_t i = 0; i < nTracks; ++i) {
	205	AliPicoTrack picotrack = static_cast<AliPicoTrack>(fTracks->At(i));
	206	if (!picotrack)
	207	continue;
	208
	209	Bool_t accept = kTRUE;
	210	Double_t eff[3] = {0};
	211	Double_t rnd = fRandom->Uniform(1.);
	212	if(fUseDiceEfficiency) accept = DiceEfficiency(picotrack,eff,rnd);
	213	if(!accept) continue;
	214
	215	AliPicoTrack *track = NULL;
	216	if(fUseTrPtResolutionSmearing) {
	217	track = SmearPt(picotrack,eff,rnd);
4ccfe90c	218	(*fTracksOut)[it] = track;
b3821b05	219	} else
4ccfe90c	220	track = new ((fTracksOut)[it]) AliPicoTrack(picotrack);
b3821b05	221
	222	track->SetBit(TObject::kBitMask,1);
	223	fHistPtDet->Fill(track->Pt());
4ccfe90c	224	it++;
b3821b05	225	}
	226	}
	227
	228	//________________________________________________________________________
	229	Bool_t AliJetFastSimulation::DiceEfficiency(AliPicoTrack *vp, Double_t eff[3], Double_t rnd)
	230	{
	231	// Dice to decide if particle is kept or not - toy model for efficiency
	232	//
	233	Double_t sumEff = 0.;
	234	Double_t pT = 0.;
	235
	236	if(fEfficiencyFixed<1.)
	237	sumEff = fEfficiencyFixed;
	238	else {
	239	pT = vp->Pt();
	240	Double_t pTtmp = pT;
	241	if(pT>10.) pTtmp = 10.;
	242	if(fhEffH1) eff[0] = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
	243	if(fhEffH2) eff[1] = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
	244	if(fhEffH3) eff[2] = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
	245
	246	sumEff = eff[0]+eff[1]+eff[2];
	247	}
	248	fp1Efficiency->Fill(vp->Pt(),sumEff);
	249	if(rnd>sumEff && pT > fDiceEfficiencyMinPt) return kFALSE;
	250	return kTRUE;
	251	}
	252
	253	//________________________________________________________________________
	254	AliPicoTrack* AliJetFastSimulation::SmearPt(AliPicoTrack *vp, Double_t eff[3], Double_t rnd)
	255	{
	256	//Smear momentum of generated particle
	257	Double_t smear = 1.;
	258	Double_t pT = vp->Pt();
	259	//Select hybrid track category
	260
	261	//Sort efficiencies from large to small
	262	Int_t cat[3] = {0};
	263	TMath::Sort(3,eff,cat);
	264	if(rnd<=eff[cat[2]])
	265	smear = GetMomentumSmearing(cat[2],pT);
	266	else if(rnd<=(eff[cat[2]]+eff[cat[1]]))
	267	smear = GetMomentumSmearing(cat[1],pT);
	268	else
	269	smear = GetMomentumSmearing(cat[0],pT);
	270
	271	fp1PtResolution->Fill(vp->Pt(),smear);
	272
	273	Double_t sigma = vp->Pt()*smear;
	274	Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
	275	fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
	276
	277	AliPicoTrack *picotrack = new AliPicoTrack(pTrec,
	278	vp->Eta(),
	279	vp->Phi(),
	280	vp->Charge(),
	281	vp->GetLabel(),
	282	AliPicoTrack::GetTrackType(vp),
	283	vp->GetTrackEtaOnEMCal(),
	284	vp->GetTrackPhiOnEMCal(),
	285	vp->GetTrackPtOnEMCal(),
	286	vp->IsEMCAL(),
	287	0.13957); //assume pion mass
	288	return picotrack;
289	}
290
291	//________________________________________________________________________
292	Double_t AliJetFastSimulation::GetMomentumSmearing(Int_t cat, Double_t pt) {
293
294	//
295	// Get smearing on generated momentum
296	//
297
298	TProfile *fMomRes = 0x0;
299	if(cat==1 && fMomResH1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
300	if(cat==2 && fMomResH2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
301	if(cat==3 && fMomResH3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
302
303	if(!fMomRes)
304	return 0.;
305
306	Double_t smear = 0.;
307	if(pt>20.) {
308	if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
309	if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
310	if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
311	}
312	else {
313	Int_t bin = fMomRes->FindBin(pt);
314	smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
315	}
316
317	if(fMomRes) delete fMomRes;
318
319	return smear;
320	}
321
322	//________________________________________________________________________
323	void AliJetFastSimulation::LoadTrPtResolutionRootFileFromOADB() {
324
325	if (!gGrid && fPathTrPtResolution.Contains("alien://")) {
326	AliInfo("Trying to connect to AliEn ...");
327	TGrid::Connect("alien://");
328	}
329
330	TFile *f = TFile::Open(fPathTrPtResolution.Data());
331	if(!f)return;
332	TProfile *fProfPtPtSigma1PtGlobSt = NULL;
333	TProfile *fProfPtPtSigma1PtGlobCnoSPD = NULL;
334	TProfile *fProfPtPtSigma1PtGlobCnoITS = NULL;
335	if(fNTrackClasses>0) fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
336	if(fNTrackClasses>1) fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
337	if(fNTrackClasses>2) fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
338
339	SetSmearResolution(kTRUE);
340	SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoSPD,fProfPtPtSigma1PtGlobCnoITS);
341	}
342
343	//________________________________________________________________________
344	void AliJetFastSimulation::LoadTrEfficiencyRootFileFromOADB() {
345
346	if (!gGrid && fPathTrPtResolution.Contains("alien://")) {
347	AliInfo("Trying to connect to AliEn ...");
348	TGrid::Connect("alien://");
349	}
350
351	TFile *f = TFile::Open(fPathTrEfficiency.Data());
352	if(!f)return;
353	TH1D *hEffPosGlobSt = NULL;
354	TH1D *hEffPosGlobCnoSPD = NULL;
355	TH1D *hEffPosGlobCnoITS = NULL;
356	if(fNTrackClasses>0) hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
357	if(fNTrackClasses>1) hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
358	if(fNTrackClasses>2) hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
359
360	SetDiceEfficiency(kTRUE);
361	SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoSPD,hEffPosGlobCnoITS);
362	}
363
364	//________________________________________________________________________
365	void AliJetFastSimulation::SetMomentumResolutionHybrid(TProfile p1, TProfile p2, TProfile *p3) {
366	//
367	// set mom res profiles
368	//
369	if(fMomResH1) delete fMomResH1;
370	if(fMomResH2) delete fMomResH2;
371	if(fMomResH3) delete fMomResH3;
372
373	if(p1) fMomResH1 = new TProfile(p1);//(TProfile)p1->Clone("fMomResH1");
374	if(p2) fMomResH2 = new TProfile(p2);//(TProfile)p2->Clone("fMomResH2");
375	if(p3) fMomResH3 = new TProfile(p3);//(TProfile)p3->Clone("fMomResH3");
376	}
377
378	//________________________________________________________________________
379	void AliJetFastSimulation:: SetEfficiencyHybrid(TH1 h1, TH1 h2, TH1 *h3) {
380	//
381	// set tracking efficiency histos
382	//
383	if(h1) fhEffH1 = (TH1*)h1->Clone("fhEffH1");
384	if(h2) fhEffH2 = (TH1*)h2->Clone("fhEffH2");
385	if(h3) fhEffH3 = (TH1*)h3->Clone("fhEffH3");
386	}
387
388	//________________________________________________________________________
389	void AliJetFastSimulation::FitMomentumResolution() {
390	//
391	// Fit linear function on momentum resolution at high pT
392	//
393
394	if(!fMomResH1Fit && fMomResH1) {
395	fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
396	fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
397	fMomResH1Fit ->SetRange(5.,100.);
398	}
399
400	if(!fMomResH2Fit && fMomResH2) {
401	fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
402	fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
403	fMomResH2Fit ->SetRange(5.,100.);
404	}
405
406	if(!fMomResH3Fit && fMomResH3) {
407	fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
408	fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
409	fMomResH3Fit ->SetRange(5.,100.);
410	}
411
412	}
413
414