PWGJE/UserTasks/AliAnalysisTaskJetHBOM.cxx

   1 // **************************************
   2 // Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
   3 // *******************************************
   4
   5
   6 /**************************************************************************
   7  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   8  *                                                                        *
   9  * Author: The ALICE Off-line Project.                                    *
  10  * Contributors are mentioned in the code where appropriate.              *
  11  *                                                                        *
  12  * Permission to use, copy, modify and distribute this software and its   *
  13  * documentation strictly for non-commercial purposes is hereby granted   *
  14  * without fee, provided that the above copyright notice appears in all   *
  15  * copies and that both the copyright notice and this permission notice   *
  16  * appear in the supporting documentation. The authors make no claims     *
  17  * about the suitability of this software for any purpose. It is          *
  18  * provided "as is" without express or implied warranty.                  *
  19  **************************************************************************/
  20
  21
  22 #include <TROOT.h>
  23 #include <TRandom3.h>
  24 #include <TSystem.h>
  25 #include <TInterpreter.h>
  26 #include <TChain.h>
  27 #include <TRefArray.h>
  28 #include <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TF1.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetHBOM.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDEvent.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODHandler.h"
  45 #include "AliAODExtension.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODJet.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliMCEvent.h"
  51 #include "AliStack.h"
  52 #include "AliGenPythiaEventHeader.h"
  53 #include "AliGenCocktailEventHeader.h"
  54 #include "AliInputEventHandler.h"
  55 #include "AliAODJetEventBackground.h"
  56
  57
  58 using std::cout;
  59 using std::endl;
  60 using std::vector;
  61
  62 ClassImp(AliAnalysisTaskJetHBOM)
  63
  64 AliAnalysisTaskJetHBOM::~AliAnalysisTaskJetHBOM(){
  65   //
  66   // Destructor
  67   //
  68
  69   delete fRef;
  70   delete fRandom;
  71
  72   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
  73   delete fTCARandomConesOut;
  74
  75   //if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
  76   //delete fTCARandomConesOutRan;
  77
  78 }
  79
  80 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM():
  81   AliAnalysisTaskSE(),
  82   fAOD(0x0),
  83   fAODExtension(0x0),
  84   fRef(new TRefArray),
  85   fUseAODTrackInput(kFALSE),
  86   fUseAODMCInput(kFALSE),
  87   fEventSelection(kFALSE),
  88   fFilterMask(0),
  89   fFilterMaskBestPt(0),
  90   fFilterType(0),
  91   fJetTypes(kJet),
  92   fTrackTypeRec(kTrackUndef),
  93   fTrackTypeGen(kTrackUndef),
  94   fNSkipLeadingRan(0),
  95   fNSkipLeadingCone(0),
  96   fNRandomCones(0),
  97   randCone_pos(0),
  98   randCone_Eta(0),
  99   randCone_Phi(0),
 100   fNHBOM(0),
 101   fTrackEtaWindow(0.9),
 102   //fRecEtaWindow(0.5),
 103   fTrackPtCut(0.),
 104   fJetOutputMinPt(0.150),
 105   fMaxTrackPtInJet(100.),
 106   fVtxZCut(8),
 107   fVtxR2Cut(1),
 108   fCentCutUp(0),
 109   fCentCutLo(0),
 110   fNonStdBranch(""),
 111   fBackgroundBranch(""),
 112   fNonStdFile(""),
 113   fMomResH1(0x0),
 114   fMomResH2(0x0),
 115   fMomResH3(0x0),
 116   fMomResH1Fit(0x0),
 117   fMomResH2Fit(0x0),
 118   fMomResH3Fit(0x0),
 119   fhEffH1(0x0),
 120   fhEffH2(0x0),
 121   fhEffH3(0x0),
 122   fUseTrMomentumSmearing(kFALSE),
 123   fUseDiceEfficiency(kFALSE),
 124   fRparam(1.0),
 125   fAlgorithm(fastjet::kt_algorithm),
 126   fStrategy(fastjet::Best),
 127   fRecombScheme(fastjet::BIpt_scheme),
 128   fAreaType(fastjet::active_area),
 129   fGhostArea(0.01),
 130   fActiveAreaRepeats(1),
 131   fGhostEtamax(1.5),
 132   background(0),
 133   fTCARandomConesOut(0x0),
 134   //fTCARandomConesOutRan(0x0),
 135   fRandom(0),
 136   fh1Xsec(0x0),
 137   fh1Trials(0x0),
 138   fh1PtHard(0x0),
 139   fh1PtHardNoW(0x0),
 140   fh1PtHardTrials(0x0),
 141   fh1Nch(0x0),
 142   fh1CentralityPhySel(0x0),
 143   fh1Centrality(0x0),
 144   fh1DeltapT(0x0),
 145   fh1Rho(0x0),
 146   fh1PtRandCone(0x0),
 147   //fh1Area(0x0),
 148   fh1efficiencyPt(0x0),
 149   fh2efficiencyPhi(0x0),
 150   fh1ZPhySel(0x0),
 151   fh1Z(0x0),
 152   fHistList(0x0)
 153 {
 154
 155 }
 156
 157 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(const char* name):
 158   AliAnalysisTaskSE(name),
 159   fAOD(0x0),
 160   fAODExtension(0x0),
 161   fRef(new TRefArray),
 162   fUseAODTrackInput(kFALSE),
 163   fUseAODMCInput(kFALSE),
 164   fEventSelection(kFALSE),
 165   fFilterMask(0),
 166   fFilterMaskBestPt(0),
 167   fFilterType(0),
 168   fJetTypes(kJet),
 169   fTrackTypeRec(kTrackUndef),
 170   fTrackTypeGen(kTrackUndef),
 171   fNSkipLeadingRan(0),
 172   fNSkipLeadingCone(0),
 173   fNRandomCones(0),
 174   randCone_pos(0),
 175   randCone_Eta(0),
 176   randCone_Phi(0),
 177   fNHBOM(0),
 178   fTrackEtaWindow(0.9),
 179   //fRecEtaWindow(0.5),
 180   fTrackPtCut(0.),
 181   fJetOutputMinPt(0.150),
 182   fMaxTrackPtInJet(100.),
 183   fVtxZCut(8),
 184   fVtxR2Cut(1),
 185   fCentCutUp(0),
 186   fCentCutLo(0),
 187   fNonStdBranch(""),
 188   fBackgroundBranch(""),
 189   fNonStdFile(""),
 190   fMomResH1(0x0),
 191   fMomResH2(0x0),
 192   fMomResH3(0x0),
 193   fMomResH1Fit(0x0),
 194   fMomResH2Fit(0x0),
 195   fMomResH3Fit(0x0),
 196   fhEffH1(0x0),
 197   fhEffH2(0x0),
 198   fhEffH3(0x0),
 199   fUseTrMomentumSmearing(kFALSE),
 200   fUseDiceEfficiency(kFALSE),
 201   fRparam(1.0),
 202   fAlgorithm(fastjet::kt_algorithm),
 203   fStrategy(fastjet::Best),
 204   fRecombScheme(fastjet::BIpt_scheme),
 205   fAreaType(fastjet::active_area),
 206   fGhostArea(0.01),
 207   fActiveAreaRepeats(1),
 208   fGhostEtamax(1.5),
 209   background(0),
 210   fTCARandomConesOut(0x0),
 211   //fTCARandomConesOutRan(0x0),
 212   fRandom(0),
 213   fh1Xsec(0x0),
 214   fh1Trials(0x0),
 215   fh1PtHard(0x0),
 216   fh1PtHardNoW(0x0),
 217   fh1PtHardTrials(0x0),
 218   fh1Nch(0x0),
 219   fh1CentralityPhySel(0x0),
 220   fh1Centrality(0x0),
 221   fh1DeltapT(0x0),
 222   fh1Rho(0x0),
 223   fh1PtRandCone(0x0),
 224   //fh1Area(0x0),
 225   fh1efficiencyPt(0x0),
 226   fh2efficiencyPhi(0x0),
 227   fh1ZPhySel(0x0),
 228   fh1Z(0x0),
 229   fHistList(0x0)
 230 {
 231  DefineOutput(1, TList::Class());
 232 }
 233
 234
 235
 236 Bool_t AliAnalysisTaskJetHBOM::Notify()
 237 {
 238   //
 239   // Implemented Notify() to read the cross sections
 240   // and number of trials from pyxsec.root
 241   //
 242   return kTRUE;
 243 }
 244
 245 void AliAnalysisTaskJetHBOM::UserCreateOutputObjects()
 246 {
 247
 248   //
 249   // Create the output container
 250   //
 251
 252   fRandom = new TRandom3(0);
 253
 254
 255   // Connect the AOD
 256
 257
 258   if (fDebug > 1) printf("AnalysisTaskJetHBOM::UserCreateOutputObjects() \n");
 259
 260
 261
 262   if(fNonStdBranch.Length()!=0)
 263     {
 264       // only create the output branch if we have a name
 265       // Create a new branch for jets...
 266       //  -> cleared in the UserExec....
 267       // here we can also have the case that the brnaches are written to a separate file
 268
 269       // create the branch for the random cones with the same R
 270       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
 271       if(fUseDiceEfficiency || fUseTrMomentumSmearing)
 272         cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
 273
 274       //create array for the random cones; Until now only one cone per event is used
 275       if(!AODEvent()->FindListObject(cName.Data())){
 276         fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
 277         fTCARandomConesOut->SetName(cName.Data());
 278         AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
 279       }
 280       /*
 281       // create the branch with the random for the random cones on the random event
 282       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 283       if(!AODEvent()->FindListObject(cName.Data())){
 284         fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
 285         fTCARandomConesOutRan->SetName(cName.Data());
 286         AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
 287       }
 288       */
 289       if(fNonStdFile.Length()!=0){
 290         //
 291         // case that we have an AOD extension we need to fetch the jets from the extended output
 292         // we identify the extension aod event by looking for the branchname
 293         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 294         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 295         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 296       }
 297     }
 298
 299   //  FitMomentumResolution();
 300
 301
 302   if(!fHistList)fHistList = new TList();
 303   fHistList->SetOwner();
 304   PostData(1, fHistList); // post data in any case once
 305
 306   Bool_t oldStatus = TH1::AddDirectoryStatus();
 307   TH1::AddDirectory(kFALSE);
 308
 309   //
 310   //  Histogram
 311
 312   const Int_t nBinPt = 100;
 313   Double_t binLimitsPt[nBinPt+1];
 314   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 315     if(iPt == 0){
 316       binLimitsPt[iPt] = 0.0;
 317     }
 318     else {// 1.0
 319       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
 320     }
 321   }
 322
 323   const Int_t nBinPhi = 90;
 324   Double_t binLimitsPhi[nBinPhi+1];
 325   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 326     if(iPhi==0){
 327       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 328     }
 329     else{
 330       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 331     }
 332   }
 333
 334
 335
 336   const Int_t nBinEta = 40;
 337   Double_t binLimitsEta[nBinEta+1];
 338   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 339     if(iEta==0){
 340       binLimitsEta[iEta] = -2.0;
 341     }
 342     else{
 343       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 344     }
 345   }
 346
 347   const int nChMax = 5000;
 348
 349   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 350   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 351
 352   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 353   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 354
 355
 356   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 357   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 358   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 359
 360   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 361
 362   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 363   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 364
 365   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 366   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 367
 368   fh1DeltapT = new TH1F("fh1DeltapT","DeltapT",100,-50,50);
 369   fh1Rho = new TH1F("fh1Rho","Rho",100,0,200);
 370   fh1PtRandCone = new TH1F("fh1PtRandCone","pt",100,0,200);
 371   //fh1Area = new TH1F("fh1Area","area",100,0,1);
 372
 373   const Int_t saveLevel = 3; // large save level more histos
 374   if(saveLevel>0){
 375     fHistList->Add(fh1Xsec);
 376     fHistList->Add(fh1Trials);
 377
 378     fHistList->Add(fh1Nch);
 379     fHistList->Add(fh1Centrality);
 380     fHistList->Add(fh1CentralityPhySel);
 381     fHistList->Add(fh1Z);
 382     fHistList->Add(fh1ZPhySel);
 383     fHistList->Add(fh1DeltapT);
 384     fHistList->Add(fh1Rho);
 385     fHistList->Add(fh1PtRandCone);
 386     //fHistList->Add(fh1Area);
 387   }
 388
 389   // =========== Switch on Sumw2 for all histos ===========
 390   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 391     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 392     if (h1){
 393       h1->Sumw2();
 394       continue;
 395     }
 396     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 397     if(hn)hn->Sumw2();
 398   }
 399   TH1::AddDirectory(oldStatus);
 400 }
 401
 402 void AliAnalysisTaskJetHBOM::Init()
 403 {
 404   //
 405   // Initialization
 406   //
 407
 408   if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
 409
 410   FitMomentumResolution();
 411
 412 }
 413
 414 void AliAnalysisTaskJetHBOM::UserExec(Option_t */*option*/)
 415 {
 416
 417   // handle and reset the output jet branch
 418
 419   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
 420   //if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
 421
 422   //
 423   // Execute analysis for current event
 424   //
 425   AliESDEvent *fESD = 0;
 426   if(fUseAODTrackInput){
 427     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 428     if(!fAOD){
 429       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 430       return;
 431     }
 432     // fetch the header
 433   }
 434   else{
 435     //  assume that the AOD is in the general output...
 436     fAOD  = AODEvent();
 437     if(!fAOD){
 438       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 439       return;
 440     }
 441     if(fDebug>0){
 442       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 443     }
 444   }
 445
 446   //Check if information is provided detector level effects
 447   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
 448   if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
 449
 450   Bool_t selectEvent =  false;
 451   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 452
 453   Float_t cent = 0;
 454   Float_t zVtx  = 0;
 455
 456   if(fAOD){
 457     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 458     TString vtxTitle(vtxAOD->GetTitle());
 459     zVtx = vtxAOD->GetZ();
 460
 461     cent = fAOD->GetHeader()->GetCentrality();
 462     if(physicsSelection){
 463       fh1CentralityPhySel->Fill(cent);
 464       fh1ZPhySel->Fill(zVtx);
 465     }
 466     // zVertex and centrality selection
 467     if(fEventSelection){
 468       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 469         Float_t yvtx = vtxAOD->GetY();
 470         Float_t xvtx = vtxAOD->GetX();
 471         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 472         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){//usual fVtxZCut=10 and fVtxR2Cut=1  // apply vertex cut later on
 473           if(physicsSelection){
 474             selectEvent = true;
 475           }
 476         }
 477       }
 478       //centrality selection
 479       if(fCentCutUp>0){
 480         if(cent<fCentCutLo||cent>fCentCutUp){
 481           selectEvent = false;
 482         }
 483       }
 484     }else{
 485       selectEvent = true;
 486     }
 487   }
 488
 489
 490   if(!selectEvent){
 491     PostData(1, fHistList);
 492     return;
 493   }
 494   fh1Centrality->Fill(cent);
 495   fh1Z->Fill(zVtx);
 496   fh1Trials->Fill("#sum{ntrials}",1);
 497
 498
 499   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 500
 501   // ==== General variables needed
 502
 503
 504
 505   // we simply fetch the tracks/mc particles as a list of AliVParticles
 506
 507   //reconstructed particles
 508   TList recParticles;
 509   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 510   Float_t nCh = recParticles.GetEntries();
 511   fh1Nch->Fill(nCh);
 512   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 513   //nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 514   //if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 515
 516
 517   //apply efficiency fNHBOM times
 518   if(fNHBOM>0){
 519     for(int particle=0;particle<recParticles.GetEntries();particle++){
 520       // hier Effizienzen laden und überprüfen ob das Teilchen nachgewiesen wird.
 521       AliVParticle *vp = (AliVParticle*)recParticles.At(particle);
 522
 523       Double_t pT = vp->Pt();
 524       Double_t phi = vp->Phi();
 525
 526       //load efficiency
 527       Double_t efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
 528       Double_t efficiencyPhi = fh2efficiencyPhi->GetBinContent(fh2efficiencyPhi->FindBin(phi,pT));
 529       if(pT>10){
 530         efficiencyPt = 0.857; //this is the result for the fit with pT>10; statistic is low for pT>10
 531         //if the efficiency is not from fastMCInput_LHC10h_110719a.root this should be calculated new
 532         if(fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))>0.9||fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))<0.8){
 533           efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
 534           Printf("%s:%d Wrong efficiency input. Check efficiency for pt>10GeV",(char*)__FILE__,__LINE__);
 535         }
 536       }
 537       Double_t eff = efficiencyPt*efficiencyPhi; //over all efficiency
 538
 539       // if ran<eff -> particle is detected eff^fNHBOM = efficiency to detect particle fNHBOM times
 540       Double_t ran = fRandom->Rndm();
 541       if(ran>TMath::Power(eff,fNHBOM)){
 542         recParticles.Remove(vp);
 543       }
 544     }
 545   }
 546
 547
 548   // find the jets....
 549
 550   vector<fastjet::PseudoJet> inputParticlesRec;
 551   vector<fastjet::PseudoJet> inputParticlesRecRan;
 552
 553   //randomize particles
 554   AliAODJet vTmpRan(1,0,0,1);
 555   for(int i = 0; i < recParticles.GetEntries(); i++){
 556     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 557
 558     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 559     // we take total momentum here
 560
 561       //Add particles to fastjet in case we are not running toy model
 562       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 563       jInp.set_user_index(i);
 564       inputParticlesRec.push_back(jInp);
 565
 566     // the randomized input changes eta and phi, but keeps the p_T
 567       Double_t pT = vp->Pt();
 568       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
 569       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 570
 571       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 572       Double_t pZ = pT/TMath::Tan(theta);
 573
 574       Double_t pX = pT * TMath::Cos(phi);
 575       Double_t pY = pT * TMath::Sin(phi);
 576       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 577       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 578
 579       jInpRan.set_user_index(i);
 580       inputParticlesRecRan.push_back(jInpRan);
 581       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
 582
 583   }// recparticles
 584
 585   if(inputParticlesRec.size()==0){
 586     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 587     PostData(1, fHistList);
 588     return;
 589   }
 590
 591   // run fast jet
 592   // employ setters for these...
 593
 594   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 595   fastjet::AreaType areaType =   fastjet::active_area;
 596   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 597   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 598   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 599
 600   //range where to compute background
 601   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 602   phiMin = 0;
 603   phiMax = 2*TMath::Pi();
 604   rapMax = fGhostEtamax - fRparam;
 605   rapMin = - fGhostEtamax + fRparam;
 606   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 607
 608
 609   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
 610   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
 611
 612
 613  // loop over all jets and fill information, first one is the leading jet
 614
 615   if(inclusiveJets.size()>0){
 616
 617     //background estimates:all bckg jets wo the 2 hardest
 618     vector<fastjet::PseudoJet> jets2=sortedJets;
 619     if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); //removes the two jets with the highest pT; +2 is correct ro remove 2 jets
 620     Double_t bkg1=0;
 621     Double_t sigma1=0.;
 622     Double_t meanarea1=0.;
 623     clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
 624     background = bkg1;//sets background variable of the task to the correct value
 625
 626
 627     // generate random cones
 628     if(fTCARandomConesOut){
 629       // create a random jet within the acceptance
 630       Double_t etaMax = fTrackEtaWindow - fRparam;//0.9 - 0.4
 631       Int_t nCone = 0;
 632       //      Int_t nConeRan = 0;
 633       Double_t pTC = 1; // small number
 634       Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 635       Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 636       // use fixed position for random Cones
 637       if(randCone_pos){
 638         etaC = randCone_Eta;
 639         phiC = randCone_Phi;
 640       }
 641       // massless jet
 642       Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 643       Double_t pZC = pTC/TMath::Tan(thetaC);
 644       Double_t pXC = pTC * TMath::Cos(phiC);
 645       Double_t pYC = pTC * TMath::Sin(phiC);
 646       Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 647       AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 648
 649       tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
 650       if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
 651       //if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
 652
 653       // loop over the reconstructed particles and add up the pT in the random cones
 654       // maybe better to loop over randomized particles not in the real jets...
 655       // but this by definition brings dow average energy in the whole  event
 656       AliAODJet vTmpRanR(1,0,0,1);
 657       for(int i = 0; i < recParticles.GetEntries(); i++){
 658         AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 659         //add up energy in cone
 660         if(fTCARandomConesOut){
 661           AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(0);
 662           if(jC&&jC->DeltaR(vp)<fRparam){
 663             if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
 664             jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
 665           }
 666         }// add up energy in cone
 667         /*
 668         // the randomized input changes eta and phi, but keeps the p_T
 669         if(fTCARandomConesOutRan){
 670           Double_t pTR = vp->Pt();
 671           Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
 672           Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
 673
 674           Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
 675           Double_t pZR = pTR/TMath::Tan(thetaR);
 676
 677           Double_t pXR = pTR * TMath::Cos(phiR);
 678           Double_t pYR = pTR * TMath::Sin(phiR);
 679           Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
 680           vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
 681
 682           AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(0);
 683           if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
 684             if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
 685             jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
 686           }
 687         }
 688         */
 689       }// loop over recparticles
 690     }  //fTCARandomConesOut
 691     Float_t jetArea = fRparam*fRparam*TMath::Pi();
 692     if(fTCARandomConesOut){
 693       // rescale the momentum vectors for the random cones
 694       AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(0);
 695       if(rC){
 696         Double_t etaC = rC->Eta();
 697         Double_t phiC = rC->Phi();
 698         // massless jet, unit vector
 699         Double_t pTC = rC->ChargedBgEnergy();
 700         if(pTC<=0)pTC = 0.001; // for almost empty events
 701         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 702         Double_t pZC = pTC/TMath::Tan(thetaC);
 703         Double_t pXC = pTC * TMath::Cos(phiC);
 704         Double_t pYC = pTC * TMath::Sin(phiC);
 705         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 706         rC->SetPxPyPzE(pXC,pYC,pZC, pC);
 707         rC->SetBgEnergy(0,0);
 708         rC->SetEffArea(jetArea,0);
 709       }
 710     }
 711     /*
 712     if(fTCARandomConesOutRan){
 713       AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(0);
 714       // same with random
 715       if(rC){
 716         Double_t etaC = rC->Eta();
 717         Double_t phiC = rC->Phi();
 718         // massless jet, unit vector
 719         Double_t pTC = rC->ChargedBgEnergy();
 720         if(pTC<=0)pTC = 0.001;// for almost empty events
 721         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 722         Double_t pZC = pTC/TMath::Tan(thetaC);
 723         Double_t pXC = pTC * TMath::Cos(phiC);
 724         Double_t pYC = pTC * TMath::Sin(phiC);
 725         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 726         rC->SetPxPyPzE(pXC,pYC,pZC, pC);
 727         rC->SetBgEnergy(0,0);
 728         rC->SetEffArea(jetArea,0);
 729       }
 730     }//find the random jets
 731     */
 732   }//inclusive Jets > 0
 733
 734  //Calculate delta pT
 735  AliAODJet *randCone = (AliAODJet*)fTCARandomConesOut->At(0);
 736  if(randCone){
 737    //background is the backbround density per area and area=pi*0.4^2 -> backgroundCone is the background energy under the cone
 738    Float_t backgroundCone = background * randCone->EffectiveAreaCharged();
 739    //calculates difference between expected and measured energy density
 740    Float_t ptSub = randCone->Pt() - backgroundCone;
 741    fh1DeltapT->Fill(ptSub);// delta pT
 742    fh1Rho->Fill(background);// background rho
 743    fh1PtRandCone->Fill(randCone->Pt());// pT of random cone
 744    //fh1Area->Fill(randCone->EffectiveAreaCharged()); // area of random cone; should always be pi*0.4^2 = 0.5
 745  }else{
 746    if(fDebug)Printf("%s:%d No random Cone found",(char*)__FILE__,__LINE__);
 747  }
 748
 749
 750  if (fDebug > 2){
 751    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
 752    //if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
 753  }
 754  PostData(1, fHistList);
 755 }
 756
 757 void AliAnalysisTaskJetHBOM::Terminate(Option_t */*option*/)
 758 {
 759   //
 760   // Terminate analysis
 761   //
 762     if (fDebug > 1) printf("AnalysisJetHBOM: Terminate() \n");
 763
 764     if(fMomResH1Fit) delete fMomResH1Fit;
 765     if(fMomResH2Fit) delete fMomResH2Fit;
 766     if(fMomResH3Fit) delete fMomResH3Fit;
 767
 768 }
 769
 770
 771 Int_t  AliAnalysisTaskJetHBOM::GetListOfTracks(TList *list,Int_t type){
 772
 773   //
 774   // get list of tracks/particles for different types
 775   //
 776
 777   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 778
 779   Int_t iCount = 0;
 780    if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
 781     if(type!=kTrackAODextraonly) {
 782       AliAODEvent *aod = 0;
 783       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 784       else aod = AODEvent();
 785       if(!aod){
 786         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 787         return iCount;
 788       }
 789
 790       for(int it = 0;it < aod->GetNumberOfTracks();++it){
 791         AliAODTrack *tr = aod->GetTrack(it);
 792         Bool_t bGood = false;
 793         if(fFilterType == 0)bGood = true;
 794         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
 795         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
 796         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
 797           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
 798           continue;
 799         }
 800         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
 801           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
 802           continue;
 803         }
 804         if(tr->Pt()<fTrackPtCut){
 805           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
 806           continue;
 807         }
 808         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
 809         list->Add(tr);
 810         iCount++;
 811       }
 812     }
 813     if(type==kTrackAODextra || type==kTrackAODextraonly) {
 814       AliAODEvent *aod = 0;
 815       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 816       else aod = AODEvent();
 817
 818       if(!aod){
 819         return iCount;
 820       }
 821       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
 822       if(!aodExtraTracks)return iCount;
 823       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
 824         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
 825         if (!track) continue;
 826
 827         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
 828         if(!trackAOD)continue;
 829         Bool_t bGood = false;
 830         if(fFilterType == 0)bGood = true;
 831         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
 832         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
 833         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
 834         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
 835         if(trackAOD->Pt()<fTrackPtCut) continue;
 836         list->Add(trackAOD);
 837         iCount++;
 838       }
 839     }
 840   }
 841   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 842     AliMCEvent* mcEvent = MCEvent();
 843     if(!mcEvent)return iCount;
 844     // we want to have alivpartilces so use get track
 845     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 846       if(!mcEvent->IsPhysicalPrimary(it))continue;
 847       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 848       if(type == kTrackKineAll){
 849         if(part->Pt()<fTrackPtCut)continue;
 850         list->Add(part);
 851         iCount++;
 852       }
 853       else if(type == kTrackKineCharged){
 854         if(part->Particle()->GetPDG()->Charge()==0)continue;
 855         if(part->Pt()<fTrackPtCut)continue;
 856         list->Add(part);
 857         iCount++;
 858       }
 859     }
 860   }
 861   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 862     AliAODEvent *aod = 0;
 863     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 864     else aod = AODEvent();
 865     if(!aod)return iCount;
 866     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 867     if(!tca)return iCount;
 868     for(int it = 0;it < tca->GetEntriesFast();++it){
 869       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
 870       if(!part->IsPhysicalPrimary())continue;
 871       if(type == kTrackAODMCAll){
 872         if(part->Pt()<fTrackPtCut)continue;
 873         list->Add(part);
 874         iCount++;
 875       }
 876       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 877         if(part->Charge()==0)continue;
 878         if(part->Pt()<fTrackPtCut)continue;
 879         if(kTrackAODMCCharged){
 880           list->Add(part);
 881         }
 882         else {
 883           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
 884           list->Add(part);
 885         }
 886         iCount++;
 887       }
 888     }
 889   }// AODMCparticle
 890   list->Sort();
 891   return iCount;
 892 }
 893
 894 void AliAnalysisTaskJetHBOM::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
 895
 896   //
 897   // set mom res profiles
 898   //
 899
 900   fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
 901   fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
 902   fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
 903 }
 904
 905 void AliAnalysisTaskJetHBOM:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
 906   //
 907   // set tracking efficiency histos
 908   //
 909
 910   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
 911   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
 912   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
 913 }
 914
 915 Double_t AliAnalysisTaskJetHBOM::GetMomentumSmearing(Int_t cat, Double_t pt) {
 916
 917   //
 918   // Get smearing on generated momentum
 919   //
 920
 921   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
 922
 923   TProfile *fMomRes = 0x0;
 924   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
 925   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
 926   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
 927
 928   if(!fMomRes) {
 929     return 0.;
 930   }
 931
 932
 933   Double_t smear = 0.;
 934
 935   if(pt>20.) {
 936     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
 937     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
 938     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
 939   }
 940   else {
 941
 942     Int_t bin = fMomRes->FindBin(pt);
 943
 944     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
 945
 946   }
 947
 948   if(fMomRes) delete fMomRes;
 949
 950   return smear;
 951 }
 952
 953 void AliAnalysisTaskJetHBOM::FitMomentumResolution() {
 954   //
 955   // Fit linear function on momentum resolution at high pT
 956   //
 957
 958   if(!fMomResH1Fit && fMomResH1) {
 959     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
 960     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
 961     fMomResH1Fit ->SetRange(5.,100.);
 962   }
 963
 964   if(!fMomResH2Fit && fMomResH2) {
 965     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
 966     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
 967     fMomResH2Fit ->SetRange(5.,100.);
 968   }
 969
 970   if(!fMomResH3Fit && fMomResH3) {
 971     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
 972     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
 973     fMomResH3Fit ->SetRange(5.,100.);
 974   }
 975
 976 }
 977