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