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