PWGJE/EMCALJetTasks/AliJetFastSimulation.cxx

   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),
  44   fEfficiencyFixed(1.1),
  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),
  77   fEfficiencyFixed(1.1),
  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
 182   if(fEfficiencyFixed < 1.)
 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
 193   fTracksOut->Delete();
 194   SimulateTracks();
 195   return kTRUE;
 196 }
 197
 198 //________________________________________________________________________
 199 void AliJetFastSimulation::SimulateTracks()
 200 {
 201   //Apply toy detector simulation to tracks
 202   Int_t it = 0;
 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);
 218       (*fTracksOut)[it] = track;
 219     } else
 220       track = new ((*fTracksOut)[it]) AliPicoTrack(*picotrack);
 221
 222     track->SetBit(TObject::kBitMask,1);
 223     fHistPtDet->Fill(track->Pt());
 224     it++;
 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