JETAN/DEV/AliAnalysisTaskJetCluster.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom3.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TRefArray.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TF1.h>
  36 #include <TList.h>
  37 #include <TLorentzVector.h>
  38 #include <TClonesArray.h>
  39 #include  "TDatabasePDG.h"
  40
  41 #include "AliAnalysisTaskJetCluster.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliJetFinder.h"
  44 #include "AliJetHeader.h"
  45 #include "AliJetReader.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODExtension.h"
  50 #include "AliAODTrack.h"
  51 #include "AliAODJet.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliJetKineReaderHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliInputEventHandler.h"
  60 #include "AliAODJetEventBackground.h"
  61
  62 #include "fastjet/PseudoJet.hh"
  63 #include "fastjet/ClusterSequenceArea.hh"
  64 #include "fastjet/AreaDefinition.hh"
  65 #include "fastjet/JetDefinition.hh"
  66 // get info on how fastjet was configured
  67 #include "fastjet/config.h"
  68
  69 using std::vector;
  70
  71 ClassImp(AliAnalysisTaskJetCluster)
  72
  73 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
  74   //
  75   // Destructor
  76   //
  77
  78   delete fRef;
  79   delete fRandom;
  80
  81   if(fTCAJetsOut)fTCAJetsOut->Delete();
  82   delete fTCAJetsOut;
  83
  84   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
  85   delete fTCAJetsOutRan;
  86
  87   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
  88   delete fTCARandomConesOut;
  89
  90   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
  91   delete fTCARandomConesOutRan;
  92
  93   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
  94   delete fAODJetBackgroundOut;
  95 }
  96
  97 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
  98   AliAnalysisTaskSE(),
  99   fAOD(0x0),
 100   fAODExtension(0x0),
 101   fRef(new TRefArray),
 102   fUseAODTrackInput(kFALSE),
 103   fUseAODMCInput(kFALSE),
 104   fUseBackgroundCalc(kFALSE),
 105   fEventSelection(kFALSE),
 106   fFilterMask(0),
 107   fFilterMaskBestPt(0),
 108   fFilterType(0),
 109   fJetTypes(kJet),
 110   fTrackTypeRec(kTrackUndef),
 111   fTrackTypeGen(kTrackUndef),
 112   fNSkipLeadingRan(0),
 113   fNSkipLeadingCone(0),
 114   fNRandomCones(0),
 115   fAvgTrials(1),
 116   fExternalWeight(1),
 117   fTrackEtaWindow(0.9),
 118   fRecEtaWindow(0.5),
 119   fTrackPtCut(0.),
 120   fJetOutputMinPt(0.150),
 121   fMaxTrackPtInJet(100.),
 122   fJetTriggerPtCut(0),
 123   fVtxZCut(8),
 124   fVtxR2Cut(1),
 125   fCentCutUp(0),
 126   fCentCutLo(0),
 127   fNonStdBranch(""),
 128   fBackgroundBranch(""),
 129   fNonStdFile(""),
 130   fMomResH1(0x0),
 131   fMomResH2(0x0),
 132   fMomResH3(0x0),
 133   fMomResH1Fit(0x0),
 134   fMomResH2Fit(0x0),
 135   fMomResH3Fit(0x0),
 136   fhEffH1(0x0),
 137   fhEffH2(0x0),
 138   fhEffH3(0x0),
 139   fUseTrPtResolutionSmearing(kFALSE),
 140   fUseDiceEfficiency(kFALSE),
 141   fDiceEfficiencyMinPt(0.),
 142   fUseTrPtResolutionFromOADB(kFALSE),
 143   fUseTrEfficiencyFromOADB(kFALSE),
 144   fPathTrPtResolution(""),
 145   fPathTrEfficiency(""),
 146   fChangeEfficiencyFraction(0.),
 147   fRparam(1.0),
 148   fAlgorithm(fastjet::kt_algorithm),
 149   fStrategy(fastjet::Best),
 150   fRecombScheme(fastjet::BIpt_scheme),
 151   fAreaType(fastjet::active_area),
 152   fGhostArea(0.01),
 153   fActiveAreaRepeats(1),
 154   fGhostEtamax(1.5),
 155   fTCAJetsOut(0x0),
 156   fTCAJetsOutRan(0x0),
 157   fTCARandomConesOut(0x0),
 158   fTCARandomConesOutRan(0x0),
 159   fAODJetBackgroundOut(0x0),
 160   fRandom(0),
 161   fh1Xsec(0x0),
 162   fh1Trials(0x0),
 163   fh1PtHard(0x0),
 164   fh1PtHardNoW(0x0),
 165   fh1PtHardTrials(0x0),
 166   fh1NJetsRec(0x0),
 167   fh1NConstRec(0x0),
 168   fh1NConstLeadingRec(0x0),
 169   fh1PtJetsRecIn(0x0),
 170   fh1PtJetsLeadingRecIn(0x0),
 171   fh1PtJetConstRec(0x0),
 172   fh1PtJetConstLeadingRec(0x0),
 173   fh1PtTracksRecIn(0x0),
 174   fh1PtTracksLeadingRecIn(0x0),
 175   fh1NJetsRecRan(0x0),
 176   fh1NConstRecRan(0x0),
 177   fh1PtJetsLeadingRecInRan(0x0),
 178   fh1NConstLeadingRecRan(0x0),
 179   fh1PtJetsRecInRan(0x0),
 180   fh1PtTracksGenIn(0x0),
 181   fh1Nch(0x0),
 182   fh1CentralityPhySel(0x0),
 183   fh1Centrality(0x0),
 184   fh1CentralitySelect(0x0),
 185   fh1ZPhySel(0x0),
 186   fh1Z(0x0),
 187   fh1ZSelect(0x0),
 188   fh2NRecJetsPt(0x0),
 189   fh2NRecTracksPt(0x0),
 190   fh2NConstPt(0x0),
 191   fh2NConstLeadingPt(0x0),
 192   fh2JetPhiEta(0x0),
 193   fh2LeadingJetPhiEta(0x0),
 194   fh2JetEtaPt(0x0),
 195   fh2LeadingJetEtaPt(0x0),
 196   fh2TrackEtaPt(0x0),
 197   fh2LeadingTrackEtaPt(0x0),
 198   fh2JetsLeadingPhiEta(0x0),
 199   fh2JetsLeadingPhiPt(0x0),
 200   fh2TracksLeadingPhiEta(0x0),
 201   fh2TracksLeadingPhiPt(0x0),
 202   fh2TracksLeadingJetPhiPt(0x0),
 203   fh2JetsLeadingPhiPtW(0x0),
 204   fh2TracksLeadingPhiPtW(0x0),
 205   fh2TracksLeadingJetPhiPtW(0x0),
 206   fh2NRecJetsPtRan(0x0),
 207   fh2NConstPtRan(0x0),
 208   fh2NConstLeadingPtRan(0x0),
 209   fh2PtNch(0x0),
 210   fh2PtNchRan(0x0),
 211   fh2PtNchN(0x0),
 212   fh2PtNchNRan(0x0),
 213   fh2TracksLeadingJetPhiPtRan(0x0),
 214   fh2TracksLeadingJetPhiPtWRan(0x0),
 215   fh3CentvsRhoLeadingTrackCorr(0x0),
 216   fh3CentvsSigmaLeadingTrackCorr(0x0),
 217   fh3MultvsRhoLeadingTrackCorr(0x0),
 218   fh3MultvsSigmaLeadingTrackCorr(0x0),
 219   fh2PtGenPtSmeared(0x0),
 220   fp1Efficiency(0x0),
 221   fp1PtResolution(0x0),
 222   fHistList(0x0)
 223 {
 224   //
 225   // Constructor
 226   //
 227
 228   for(int i = 0;i<3;i++){
 229     fh1BiARandomCones[i] = 0;
 230     fh1BiARandomConesRan[i] = 0;
 231   }
 232   for(int i = 0;i<kMaxCent;i++){
 233     fh2JetsLeadingPhiPtC[i] = 0;
 234     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 235     fh2TracksLeadingJetPhiPtC[i] = 0;
 236     fh2TracksLeadingJetPhiPtWC[i] = 0;
 237   }
 238 }
 239
 240 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 241   AliAnalysisTaskSE(name),
 242   fAOD(0x0),
 243   fAODExtension(0x0),
 244   fRef(new TRefArray),
 245   fUseAODTrackInput(kFALSE),
 246   fUseAODMCInput(kFALSE),
 247   fUseBackgroundCalc(kFALSE),
 248   fEventSelection(kFALSE),                                                        fFilterMask(0),
 249   fFilterMaskBestPt(0),
 250   fFilterType(0),
 251   fJetTypes(kJet),
 252   fTrackTypeRec(kTrackUndef),
 253   fTrackTypeGen(kTrackUndef),
 254   fNSkipLeadingRan(0),
 255   fNSkipLeadingCone(0),
 256   fNRandomCones(0),
 257   fAvgTrials(1),
 258   fExternalWeight(1),
 259   fTrackEtaWindow(0.9),
 260   fRecEtaWindow(0.5),
 261   fTrackPtCut(0.),
 262   fJetOutputMinPt(0.150),
 263   fMaxTrackPtInJet(100.),
 264   fJetTriggerPtCut(0),
 265   fVtxZCut(8),
 266   fVtxR2Cut(1),
 267   fCentCutUp(0),
 268   fCentCutLo(0),
 269   fNonStdBranch(""),
 270   fBackgroundBranch(""),
 271   fNonStdFile(""),
 272   fMomResH1(0x0),
 273   fMomResH2(0x0),
 274   fMomResH3(0x0),
 275   fMomResH1Fit(0x0),
 276   fMomResH2Fit(0x0),
 277   fMomResH3Fit(0x0),
 278   fhEffH1(0x0),
 279   fhEffH2(0x0),
 280   fhEffH3(0x0),
 281   fUseTrPtResolutionSmearing(kFALSE),
 282   fUseDiceEfficiency(kFALSE),
 283   fDiceEfficiencyMinPt(0.),
 284   fUseTrPtResolutionFromOADB(kFALSE),
 285   fUseTrEfficiencyFromOADB(kFALSE),
 286   fPathTrPtResolution(""),
 287   fPathTrEfficiency(""),
 288   fChangeEfficiencyFraction(0.),
 289   fRparam(1.0),
 290   fAlgorithm(fastjet::kt_algorithm),
 291   fStrategy(fastjet::Best),
 292   fRecombScheme(fastjet::BIpt_scheme),
 293   fAreaType(fastjet::active_area),
 294   fGhostArea(0.01),
 295   fActiveAreaRepeats(1),
 296   fGhostEtamax(1.5),
 297   fTCAJetsOut(0x0),
 298   fTCAJetsOutRan(0x0),
 299   fTCARandomConesOut(0x0),
 300   fTCARandomConesOutRan(0x0),
 301   fAODJetBackgroundOut(0x0),
 302   fRandom(0),
 303   fh1Xsec(0x0),
 304   fh1Trials(0x0),
 305   fh1PtHard(0x0),
 306   fh1PtHardNoW(0x0),
 307   fh1PtHardTrials(0x0),
 308   fh1NJetsRec(0x0),
 309   fh1NConstRec(0x0),
 310   fh1NConstLeadingRec(0x0),
 311   fh1PtJetsRecIn(0x0),
 312   fh1PtJetsLeadingRecIn(0x0),
 313   fh1PtJetConstRec(0x0),
 314   fh1PtJetConstLeadingRec(0x0),
 315   fh1PtTracksRecIn(0x0),
 316   fh1PtTracksLeadingRecIn(0x0),
 317   fh1NJetsRecRan(0x0),
 318   fh1NConstRecRan(0x0),
 319   fh1PtJetsLeadingRecInRan(0x0),
 320   fh1NConstLeadingRecRan(0x0),
 321   fh1PtJetsRecInRan(0x0),
 322   fh1PtTracksGenIn(0x0),
 323   fh1Nch(0x0),
 324   fh1CentralityPhySel(0x0),
 325   fh1Centrality(0x0),
 326   fh1CentralitySelect(0x0),
 327   fh1ZPhySel(0x0),
 328   fh1Z(0x0),
 329   fh1ZSelect(0x0),
 330   fh2NRecJetsPt(0x0),
 331   fh2NRecTracksPt(0x0),
 332   fh2NConstPt(0x0),
 333   fh2NConstLeadingPt(0x0),
 334   fh2JetPhiEta(0x0),
 335   fh2LeadingJetPhiEta(0x0),
 336   fh2JetEtaPt(0x0),
 337   fh2LeadingJetEtaPt(0x0),
 338   fh2TrackEtaPt(0x0),
 339   fh2LeadingTrackEtaPt(0x0),
 340   fh2JetsLeadingPhiEta(0x0),
 341   fh2JetsLeadingPhiPt(0x0),
 342   fh2TracksLeadingPhiEta(0x0),
 343   fh2TracksLeadingPhiPt(0x0),
 344   fh2TracksLeadingJetPhiPt(0x0),
 345   fh2JetsLeadingPhiPtW(0x0),
 346   fh2TracksLeadingPhiPtW(0x0),
 347   fh2TracksLeadingJetPhiPtW(0x0),
 348   fh2NRecJetsPtRan(0x0),
 349   fh2NConstPtRan(0x0),
 350   fh2NConstLeadingPtRan(0x0),
 351   fh2PtNch(0x0),
 352   fh2PtNchRan(0x0),
 353   fh2PtNchN(0x0),
 354   fh2PtNchNRan(0x0),
 355   fh2TracksLeadingJetPhiPtRan(0x0),
 356   fh2TracksLeadingJetPhiPtWRan(0x0),
 357   fh3CentvsRhoLeadingTrackCorr(0x0),
 358   fh3CentvsSigmaLeadingTrackCorr(0x0),
 359   fh3MultvsRhoLeadingTrackCorr(0x0),
 360   fh3MultvsSigmaLeadingTrackCorr(0x0),
 361   fh2PtGenPtSmeared(0x0),
 362   fp1Efficiency(0x0),
 363   fp1PtResolution(0x0),
 364   fHistList(0x0)
 365 {
 366   //
 367   // named ctor
 368   //
 369
 370   for(int i = 0;i<3;i++){
 371     fh1BiARandomCones[i] = 0;
 372     fh1BiARandomConesRan[i] = 0;
 373   }
 374   for(int i = 0;i<kMaxCent;i++){
 375     fh2JetsLeadingPhiPtC[i] = 0;
 376     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 377     fh2TracksLeadingJetPhiPtC[i] = 0;
 378     fh2TracksLeadingJetPhiPtWC[i] = 0;
 379   }
 380  DefineOutput(1, TList::Class());
 381 }
 382
 383
 384
 385 Bool_t AliAnalysisTaskJetCluster::Notify()
 386 {
 387   //
 388   // Implemented Notify() to read the cross sections
 389   // and number of trials from pyxsec.root
 390   //
 391   return kTRUE;
 392 }
 393
 394 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 395 {
 396
 397   //
 398   // Create the output container
 399   //
 400
 401   fRandom = new TRandom3(0);
 402
 403
 404   // Connect the AOD
 405
 406
 407   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 408
 409
 410
 411   if(fNonStdBranch.Length()!=0)
 412     {
 413       // only create the output branch if we have a name
 414       // Create a new branch for jets...
 415       //  -> cleared in the UserExec....
 416       // here we can also have the case that the brnaches are written to a separate file
 417
 418       if(fJetTypes&kJet){
 419         fTCAJetsOut = new TClonesArray("AliAODJet", 0);
 420         fTCAJetsOut->SetName(fNonStdBranch.Data());
 421         AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
 422       }
 423
 424       if(fJetTypes&kJetRan){
 425         fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
 426         fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 427         if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
 428           fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
 429         AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
 430       }
 431
 432       if(fUseBackgroundCalc){
 433         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 434           fAODJetBackgroundOut = new AliAODJetEventBackground();
 435           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 436           if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
 437             fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
 438
 439           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
 440         }
 441       }
 442       // create the branch for the random cones with the same R
 443       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
 444       if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
 445         cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
 446
 447       if(fNRandomCones>0){
 448         if(fJetTypes&kRC){
 449           if(!AODEvent()->FindListObject(cName.Data())){
 450             fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
 451             fTCARandomConesOut->SetName(cName.Data());
 452             AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
 453           }
 454         }
 455         // create the branch with the random for the random cones on the random event
 456         if(fJetTypes&kRCRan){
 457           cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 458           if(!AODEvent()->FindListObject(cName.Data())){
 459             fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
 460             fTCARandomConesOutRan->SetName(cName.Data());
 461             AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
 462           }
 463         }
 464       }
 465
 466       if(fNonStdFile.Length()!=0){
 467         //
 468         // case that we have an AOD extension we need to fetch the jets from the extended output
 469         // we identify the extension aod event by looking for the branchname
 470         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 471         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 472         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 473       }
 474     }
 475
 476
 477   if(!fHistList)fHistList = new TList();
 478   fHistList->SetOwner();
 479   PostData(1, fHistList); // post data in any case once
 480
 481   Bool_t oldStatus = TH1::AddDirectoryStatus();
 482   TH1::AddDirectory(kFALSE);
 483
 484   //
 485   //  Histogram
 486
 487   const Int_t nBinPt = 100;
 488   Double_t binLimitsPt[nBinPt+1];
 489   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 490     if(iPt == 0){
 491       binLimitsPt[iPt] = 0.0;
 492     }
 493     else {// 1.0
 494       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
 495     }
 496   }
 497
 498   const Int_t nBinPhi = 90;
 499   Double_t binLimitsPhi[nBinPhi+1];
 500   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 501     if(iPhi==0){
 502       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 503     }
 504     else{
 505       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 506     }
 507   }
 508
 509
 510
 511   const Int_t nBinEta = 40;
 512   Double_t binLimitsEta[nBinEta+1];
 513   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 514     if(iEta==0){
 515       binLimitsEta[iEta] = -2.0;
 516     }
 517     else{
 518       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 519     }
 520   }
 521
 522   const int nChMax = 5000;
 523
 524   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 525   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 526
 527   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 528   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 529
 530
 531   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 532   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 533
 534   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 535   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 536   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 537   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 538
 539
 540   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 541   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 542   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 543
 544   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 545   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 546   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 547   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 548   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 549   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 550   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 551   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 552   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 553   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 554
 555   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 556   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
 557   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 558
 559   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 560   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
 561   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 562
 563   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 564   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 565   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 566   //
 567
 568
 569   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 570   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 571   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 572   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 573
 574   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 575   fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 576   fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 577   fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 578
 579
 580
 581   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 582                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 583   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 584                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 585
 586   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 587                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 588   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 589                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 590
 591   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 592                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 593   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 594                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 595
 596
 597
 598   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 599                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 600   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 601                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 602   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 603                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 604   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 605                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 606   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 607                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 608   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 609                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 610
 611   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 612                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 613   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 614                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 615
 616   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 617                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 618   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 619                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 620
 621   fh3CentvsRhoLeadingTrackCorr = new TH3F("fh3CentvsRhoLeadingTrackCorr","centrality vs background density; centrality; #rho; Quadrant", 100,0.,100.,600,0.,300.,4,0.,4.);
 622   fh3CentvsSigmaLeadingTrackCorr = new TH3F("fh3CentvsSigmaLeadingTrackCorr","centrality vs backgroun sigma; centrality; #sigma(#rho); Quadrant",100,0.,100.,500,0.,50.,4,0.,4.);
 623   fHistList->Add(fh3CentvsRhoLeadingTrackCorr);
 624   fHistList->Add(fh3CentvsSigmaLeadingTrackCorr);
 625
 626   fh3MultvsRhoLeadingTrackCorr = new TH3F("fh3MultvsRhoLeadingTrackCorr","mult vs background density; N_{tracks}; #rho; Quadrant", 100,0.,5000.,600,0.,300.,4,0.,4.);
 627   fh3MultvsSigmaLeadingTrackCorr = new TH3F("fh3MultvsSigmaLeadingTrackCorr","mult vs background sigma; N_{tracks}; #sigma(#rho); Quadrant",100,0.,5000.,500,0.,50.,4,0.,4.);
 628   fHistList->Add(fh3MultvsRhoLeadingTrackCorr);
 629   fHistList->Add(fh3MultvsSigmaLeadingTrackCorr);
 630
 631   //Detector level effects histos
 632   fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 633
 634   fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
 635   fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
 636
 637   fHistList->Add(fh2PtGenPtSmeared);
 638   fHistList->Add(fp1Efficiency);
 639   fHistList->Add(fp1PtResolution);
 640
 641   if(fNRandomCones>0&&fUseBackgroundCalc){
 642     for(int i = 0;i<3;i++){
 643       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 644       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 645     }
 646   }
 647
 648   for(int i = 0;i < kMaxCent;i++){
 649     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
 650     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
 651     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
 652     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
 653   }
 654
 655   const Int_t saveLevel = 3; // large save level more histos
 656   if(saveLevel>0){
 657     fHistList->Add(fh1Xsec);
 658     fHistList->Add(fh1Trials);
 659
 660     fHistList->Add(fh1NJetsRec);
 661     fHistList->Add(fh1NConstRec);
 662     fHistList->Add(fh1NConstLeadingRec);
 663     fHistList->Add(fh1PtJetsRecIn);
 664     fHistList->Add(fh1PtJetsLeadingRecIn);
 665     fHistList->Add(fh1PtTracksRecIn);
 666     fHistList->Add(fh1PtTracksLeadingRecIn);
 667     fHistList->Add(fh1PtJetConstRec);
 668     fHistList->Add(fh1PtJetConstLeadingRec);
 669     fHistList->Add(fh1NJetsRecRan);
 670     fHistList->Add(fh1NConstRecRan);
 671     fHistList->Add(fh1PtJetsLeadingRecInRan);
 672     fHistList->Add(fh1NConstLeadingRecRan);
 673     fHistList->Add(fh1PtJetsRecInRan);
 674     fHistList->Add(fh1Nch);
 675     fHistList->Add(fh1Centrality);
 676     fHistList->Add(fh1CentralitySelect);
 677     fHistList->Add(fh1CentralityPhySel);
 678     fHistList->Add(fh1Z);
 679     fHistList->Add(fh1ZSelect);
 680     fHistList->Add(fh1ZPhySel);
 681     if(fNRandomCones>0&&fUseBackgroundCalc){
 682       for(int i = 0;i<3;i++){
 683         fHistList->Add(fh1BiARandomCones[i]);
 684         fHistList->Add(fh1BiARandomConesRan[i]);
 685       }
 686     }
 687   for(int i = 0;i < kMaxCent;i++){
 688     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
 689     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
 690     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
 691     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
 692   }
 693
 694     fHistList->Add(fh2NRecJetsPt);
 695     fHistList->Add(fh2NRecTracksPt);
 696     fHistList->Add(fh2NConstPt);
 697     fHistList->Add(fh2NConstLeadingPt);
 698     fHistList->Add(fh2PtNch);
 699     fHistList->Add(fh2PtNchRan);
 700     fHistList->Add(fh2PtNchN);
 701     fHistList->Add(fh2PtNchNRan);
 702     fHistList->Add(fh2JetPhiEta);
 703     fHistList->Add(fh2LeadingJetPhiEta);
 704     fHistList->Add(fh2JetEtaPt);
 705     fHistList->Add(fh2LeadingJetEtaPt);
 706     fHistList->Add(fh2TrackEtaPt);
 707     fHistList->Add(fh2LeadingTrackEtaPt);
 708     fHistList->Add(fh2JetsLeadingPhiEta );
 709     fHistList->Add(fh2JetsLeadingPhiPt);
 710     fHistList->Add(fh2TracksLeadingPhiEta);
 711     fHistList->Add(fh2TracksLeadingPhiPt);
 712     fHistList->Add(fh2TracksLeadingJetPhiPt);
 713     fHistList->Add(fh2JetsLeadingPhiPtW);
 714     fHistList->Add(fh2TracksLeadingPhiPtW);
 715     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 716     fHistList->Add(fh2NRecJetsPtRan);
 717     fHistList->Add(fh2NConstPtRan);
 718     fHistList->Add(fh2NConstLeadingPtRan);
 719     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 720     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 721   }
 722
 723   // =========== Switch on Sumw2 for all histos ===========
 724   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 725     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 726     if (h1){
 727       h1->Sumw2();
 728       continue;
 729     }
 730     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 731     if(hn)hn->Sumw2();
 732   }
 733   TH1::AddDirectory(oldStatus);
 734 }
 735
 736 void AliAnalysisTaskJetCluster::LocalInit()
 737 {
 738   //
 739   // Initialization
 740   //
 741
 742   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 743
 744   if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
 745   if(fUseTrEfficiencyFromOADB)   LoadTrEfficiencyRootFileFromOADB();
 746
 747
 748   FitMomentumResolution();
 749
 750 }
 751
 752 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 753 {
 754
 755   // handle and reset the output jet branch
 756
 757   if(fTCAJetsOut)fTCAJetsOut->Delete();
 758   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
 759   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
 760   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
 761   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
 762
 763   AliAODJetEventBackground* externalBackground = 0;
 764   if(!externalBackground&&fBackgroundBranch.Length()){
 765     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
 766     if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 767     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 768   }
 769   //
 770   // Execute analysis for current event
 771   //
 772   AliESDEvent *fESD = 0;
 773   if(fUseAODTrackInput){
 774     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 775     if(!fAOD){
 776       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 777       return;
 778     }
 779     // fetch the header
 780   }
 781   else{
 782     //  assume that the AOD is in the general output...
 783     fAOD  = AODEvent();
 784     if(!fAOD){
 785       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 786       return;
 787     }
 788     if(fDebug>0){
 789       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 790     }
 791   }
 792
 793   //Check if information is provided detector level effects
 794   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
 795   if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
 796
 797   Bool_t selectEvent =  false;
 798   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 799
 800   Float_t cent = 0;
 801   Float_t zVtx  = 0;
 802   Int_t cenClass = -1;
 803   if(fAOD){
 804     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 805     TString vtxTitle(vtxAOD->GetTitle());
 806     zVtx = vtxAOD->GetZ();
 807
 808     cent = fAOD->GetHeader()->GetCentrality();
 809     if(cent<10)cenClass = 0;
 810     else if(cent<30)cenClass = 1;
 811     else if(cent<50)cenClass = 2;
 812     else if(cent<80)cenClass = 3;
 813     if(physicsSelection){
 814       fh1CentralityPhySel->Fill(cent);
 815       fh1ZPhySel->Fill(zVtx);
 816     }
 817
 818     if(fEventSelection){
 819       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 820         Float_t yvtx = vtxAOD->GetY();
 821         Float_t xvtx = vtxAOD->GetX();
 822         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 823         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
 824           if(physicsSelection){
 825             selectEvent = true;
 826           }
 827         }
 828       }
 829       if(fCentCutUp>0){
 830         if(cent<fCentCutLo||cent>fCentCutUp){
 831           selectEvent = false;
 832         }
 833       }
 834     }else{
 835       selectEvent = true;
 836     }
 837   }
 838
 839
 840   if(!selectEvent){
 841     PostData(1, fHistList);
 842     return;
 843   }
 844   fh1Centrality->Fill(cent);
 845   fh1Z->Fill(zVtx);
 846   fh1Trials->Fill("#sum{ntrials}",1);
 847
 848
 849   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 850
 851   // ==== General variables needed
 852
 853
 854
 855   // we simply fetch the tracks/mc particles as a list of AliVParticles
 856
 857   TList recParticles;
 858   TList genParticles;
 859
 860   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 861   Float_t nCh = recParticles.GetEntries();
 862   fh1Nch->Fill(nCh);
 863   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 864   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 865   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 866
 867   // find the jets....
 868
 869   vector<fastjet::PseudoJet> inputParticlesRec;
 870   vector<fastjet::PseudoJet> inputParticlesRecRan;
 871
 872   // Generate the random cones
 873
 874   AliAODJet vTmpRan(1,0,0,1);
 875   for(int i = 0; i < recParticles.GetEntries(); i++){
 876     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 877
 878     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 879     // we take total momentum here
 880
 881     if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
 882       //Add particles to fastjet in case we are not running toy model
 883       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 884       jInp.set_user_index(i);
 885       inputParticlesRec.push_back(jInp);
 886     }
 887     else if(fUseDiceEfficiency) {
 888
 889       // Dice to decide if particle is kept or not - toy  model for efficiency
 890       //
 891       Double_t rnd = fRandom->Uniform(1.);
 892       Double_t pT = vp->Pt();
 893       Double_t eff[3] = {0.};
 894       Double_t pTtmp = pT;
 895       if(pT>10.) pTtmp = 10.;
 896       Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
 897       Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
 898       Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
 899       Int_t cat[3] = {0};
 900       //Sort efficiencies from large to small
 901       if(eff1>eff2 && eff1>eff3) {
 902         eff[0] = eff1;
 903         cat[0] = 1;
 904         if(eff2>eff3) {
 905           eff[1] = eff2;
 906           eff[2] = eff3;
 907           cat[1] = 2;
 908           cat[2] = 3;
 909         } else {
 910           eff[1] = eff3;
 911           eff[2] = eff2;
 912           cat[1] = 3;
 913           cat[2] = 2;
 914         }
 915       }
 916       else if(eff2>eff1 && eff2>eff3) {
 917         eff[0] = eff2;
 918         cat[0] = 2;
 919         if(eff1>eff3) {
 920           eff[1] = eff1;
 921           eff[2] = eff3;
 922           cat[1] = 1;
 923           cat[2] = 3;
 924         } else {
 925           eff[1] = eff3;
 926           eff[2] = eff1;
 927           cat[1] = 3;
 928           cat[2] = 1;
 929         }
 930       }
 931       else if(eff3>eff1 && eff3>eff2) {
 932         eff[0] = eff3;
 933         cat[0] = 3;
 934         if(eff1>eff2) {
 935           eff[1] = eff1;
 936           eff[2] = eff2;
 937           cat[1] = 1;
 938           cat[2] = 2;
 939         } else {
 940           eff[1] = eff2;
 941           eff[2] = eff1;
 942           cat[1] = 2;
 943           cat[2] = 1;
 944         }
 945       }
 946
 947       Double_t sumEff = eff[0]+eff[1]+eff[2];
 948       fp1Efficiency->Fill(vp->Pt(),sumEff);
 949       if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
 950
 951       if(fUseTrPtResolutionSmearing) {
 952         //Smear momentum of generated particle
 953         Double_t smear = 1.;
 954         //Select hybrid track category
 955         if(rnd<=eff[2])
 956           smear = GetMomentumSmearing(cat[2],pT);
 957         else if(rnd<=(eff[2]+eff[1]))
 958           smear = GetMomentumSmearing(cat[1],pT);
 959         else
 960           smear = GetMomentumSmearing(cat[0],pT);
 961
 962         fp1PtResolution->Fill(vp->Pt(),smear);
 963
 964         Double_t sigma = vp->Pt()*smear;
 965         Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
 966
 967         Double_t phi   = vp->Phi();
 968         Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
 969         Double_t pX    = pTrec * TMath::Cos(phi);
 970         Double_t pY    = pTrec * TMath::Sin(phi);
 971         Double_t pZ    = pTrec/TMath::Tan(theta);
 972         Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
 973
 974         fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
 975
 976         fastjet::PseudoJet jInp(pX,pY,pZ,p);
 977         jInp.set_user_index(i);
 978         inputParticlesRec.push_back(jInp);
 979
 980       }
 981       else {
 982         fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 983         jInp.set_user_index(i);
 984         inputParticlesRec.push_back(jInp);
 985
 986       }
 987
 988     }
 989
 990     // the randomized input changes eta and phi, but keeps the p_T
 991     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 992       Double_t pT = vp->Pt();
 993       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
 994       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 995
 996       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 997       Double_t pZ = pT/TMath::Tan(theta);
 998
 999       Double_t pX = pT * TMath::Cos(phi);
1000       Double_t pY = pT * TMath::Sin(phi);
1001       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
1002       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1003
1004       jInpRan.set_user_index(i);
1005       inputParticlesRecRan.push_back(jInpRan);
1006       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1007     }
1008
1009     // fill the tref array, only needed when we write out jets
1010     if(fTCAJetsOut){
1011       if(i == 0){
1012         fRef->Delete(); // make sure to delete before placement new...
1013         if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1014           new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1015         }
1016       }
1017       if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp);  //TRefArray does not work with toy model ...
1018     }
1019   }// recparticles
1020
1021   if(inputParticlesRec.size()==0){
1022     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1023     PostData(1, fHistList);
1024     return;
1025   }
1026
1027   // run fast jet
1028   // employ setters for these...
1029
1030
1031   // now create the object that holds info about ghosts
1032   /*
1033   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1034     // reduce CPU time...
1035     fGhostArea = 0.5;
1036     fActiveAreaRepeats = 0;
1037   }
1038   */
1039
1040   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1041   fastjet::AreaType areaType =   fastjet::active_area;
1042   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1043   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1044   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1045
1046   //range where to compute background
1047   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1048   phiMin = 0;
1049   phiMax = 2*TMath::Pi();
1050   rapMax = fGhostEtamax - fRparam;
1051   rapMin = - fGhostEtamax + fRparam;
1052   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1053
1054
1055   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1056   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1057
1058
1059   fh1NJetsRec->Fill(sortedJets.size());
1060
1061  // loop over all jets an fill information, first on is the leading jet
1062
1063   Int_t nRecOver = inclusiveJets.size();
1064   Int_t nRec     = inclusiveJets.size();
1065   if(inclusiveJets.size()>0){
1066     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1067     Double_t area = clustSeq.area(sortedJets[0]);
1068     leadingJet.SetEffArea(area,0);
1069     Float_t pt = leadingJet.Pt();
1070     Int_t nAodOutJets = 0;
1071     Int_t nAodOutTracks = 0;
1072     AliAODJet *aodOutJet = 0;
1073
1074     Int_t iCount = 0;
1075     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1076       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1077       while(pt<ptCut&&iCount<nRec){
1078         nRecOver--;
1079         iCount++;
1080         if(iCount<nRec){
1081           pt = sortedJets[iCount].perp();
1082         }
1083       }
1084       if(nRecOver<=0)break;
1085       fh2NRecJetsPt->Fill(ptCut,nRecOver);
1086     }
1087     Float_t phi = leadingJet.Phi();
1088     if(phi<0)phi+=TMath::Pi()*2.;
1089     Float_t eta = leadingJet.Eta();
1090     Float_t pTback = 0;
1091     if(externalBackground){
1092       // carefull has to be filled in a task before
1093       // todo, ReArrange to the botom
1094      pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1095     }
1096     pt = leadingJet.Pt() - pTback;
1097     // correlation of leading jet with tracks
1098     TIterator *recIter = recParticles.MakeIterator();
1099     recIter->Reset();
1100     AliVParticle *tmpRecTrack = 0;
1101     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1102       Float_t tmpPt = tmpRecTrack->Pt();
1103       // correlation
1104       Float_t tmpPhi =  tmpRecTrack->Phi();
1105       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1106       Float_t dPhi = phi - tmpPhi;
1107       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1108       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1109       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1110       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1111       if(cenClass>=0){
1112         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1113         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1114       }
1115
1116     }
1117
1118
1119     TLorentzVector vecareab;
1120     for(int j = 0; j < nRec;j++){
1121       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1122       aodOutJet = 0;
1123       nAodOutTracks = 0;
1124       Float_t tmpPt = tmpRec.Pt();
1125
1126       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1127         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1128         aodOutJet->GetRefTracks()->Clear();
1129         Double_t area1 = clustSeq.area(sortedJets[j]);
1130         aodOutJet->SetEffArea(area1,0);
1131         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1132         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1133         aodOutJet->SetVectorAreaCharged(&vecareab);
1134       }
1135
1136
1137       Float_t tmpPtBack = 0;
1138       if(externalBackground){
1139         // carefull has to be filled in a task before
1140        // todo, ReArrange to the botom
1141         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1142       }
1143       tmpPt = tmpPt - tmpPtBack;
1144       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1145
1146       fh1PtJetsRecIn->Fill(tmpPt);
1147       // Fill Spectra with constituentsemacs
1148       const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1149
1150       fh1NConstRec->Fill(constituents.size());
1151       fh2PtNch->Fill(nCh,tmpPt);
1152       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1153       fh2NConstPt->Fill(tmpPt,constituents.size());
1154       // loop over constiutents and fill spectrum
1155
1156       AliVParticle *partLead = 0;
1157       Float_t ptLead = -1;
1158
1159       for(unsigned int ic = 0; ic < constituents.size();ic++){
1160         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1161         if(!part) continue;
1162         fh1PtJetConstRec->Fill(part->Pt());
1163         if(aodOutJet){
1164           if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1165           if(part->Pt()>fMaxTrackPtInJet){
1166             aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1167           }
1168         }
1169         if(part->Pt()>ptLead){
1170           ptLead = part->Pt();
1171           partLead = part;
1172         }
1173         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1174       }
1175
1176       AliAODTrack *aodT = 0;
1177       if(partLead){
1178         if(aodOutJet){
1179           //set pT of leading constituent of jet
1180           aodOutJet->SetPtLeading(partLead->Pt());
1181           aodT = dynamic_cast<AliAODTrack*>(partLead);
1182           if(aodT){
1183             if(aodT->TestFilterBit(fFilterMaskBestPt)){
1184               aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1185             }
1186           }
1187         }
1188       }
1189
1190      // correlation
1191      Float_t tmpPhi =  tmpRec.Phi();
1192      Float_t tmpEta =  tmpRec.Eta();
1193      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1194      if(j==0){
1195        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1196        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1197        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1198        fh1NConstLeadingRec->Fill(constituents.size());
1199        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1200        continue;
1201      }
1202      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1203      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1204      Float_t dPhi = phi - tmpPhi;
1205      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1206      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1207      Float_t dEta = eta - tmpRec.Eta();
1208      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1209      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1210      if(cenClass>=0){
1211        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1212        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1213      }
1214      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1215     }// loop over reconstructed jets
1216    delete recIter;
1217
1218
1219
1220    // Add the random cones...
1221    if(fNRandomCones>0&&fTCARandomConesOut){
1222      // create a random jet within the acceptance
1223      Double_t etaMax = fTrackEtaWindow - fRparam;
1224      Int_t nCone = 0;
1225      Int_t nConeRan = 0;
1226      Double_t pTC = 1; // small number
1227      for(int ir = 0;ir < fNRandomCones;ir++){
1228        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1229        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1230        // massless jet
1231        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1232        Double_t pZC = pTC/TMath::Tan(thetaC);
1233        Double_t pXC = pTC * TMath::Cos(phiC);
1234        Double_t pYC = pTC * TMath::Sin(phiC);
1235        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1236        AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1237        bool skip = false;
1238        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1239          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1240          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1241            skip = true;
1242            break;
1243          }
1244        }
1245        // test for overlap with previous cones to avoid double counting
1246        for(int iic = 0;iic<ir;iic++){
1247          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1248          if(iicone){
1249            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1250              skip = true;
1251              break;
1252            }
1253          }
1254        }
1255        if(skip)continue;
1256        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1257        tmpRecC.SetPtLeading(-1.);
1258        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1259        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1260      }// loop over random cones creation
1261
1262
1263      // loop over the reconstructed particles and add up the pT in the random cones
1264      // maybe better to loop over randomized particles not in the real jets...
1265      // but this by definition brings dow average energy in the whole  event
1266      AliAODJet vTmpRanR(1,0,0,1);
1267      for(int i = 0; i < recParticles.GetEntries(); i++){
1268        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1269        if(fTCARandomConesOut){
1270          for(int ir = 0;ir < fNRandomCones;ir++){
1271            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1272            if(jC&&jC->DeltaR(vp)<fRparam){
1273              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1274              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1275              if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1276            }
1277          }
1278        }// add up energy in cone
1279
1280        // the randomized input changes eta and phi, but keeps the p_T
1281        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1282          Double_t pTR = vp->Pt();
1283          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1284          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1285
1286          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1287          Double_t pZR = pTR/TMath::Tan(thetaR);
1288
1289          Double_t pXR = pTR * TMath::Cos(phiR);
1290          Double_t pYR = pTR * TMath::Sin(phiR);
1291          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
1292          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1293          if(fTCARandomConesOutRan){
1294            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1295              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1296              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1297                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1298                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1299                if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1300              }
1301            }
1302          }
1303        }
1304      }// loop over recparticles
1305
1306      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1307      if(fTCARandomConesOut){
1308        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1309          // rescale the momentum vectors for the random cones
1310
1311          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1312          if(rC){
1313            Double_t etaC = rC->Eta();
1314            Double_t phiC = rC->Phi();
1315            // massless jet, unit vector
1316            pTC = rC->ChargedBgEnergy();
1317            if(pTC<=0)pTC = 0.001; // for almost empty events
1318            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1319            Double_t pZC = pTC/TMath::Tan(thetaC);
1320            Double_t pXC = pTC * TMath::Cos(phiC);
1321            Double_t pYC = pTC * TMath::Sin(phiC);
1322            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1323            rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1324            rC->SetBgEnergy(0,0);
1325            rC->SetEffArea(jetArea,0);
1326          }
1327        }
1328      }
1329      if(fTCARandomConesOutRan){
1330        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1331          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1332          // same wit random
1333          if(rC){
1334            Double_t etaC = rC->Eta();
1335            Double_t phiC = rC->Phi();
1336            // massless jet, unit vector
1337            pTC = rC->ChargedBgEnergy();
1338            if(pTC<=0)pTC = 0.001;// for almost empty events
1339            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1340            Double_t pZC = pTC/TMath::Tan(thetaC);
1341            Double_t pXC = pTC * TMath::Cos(phiC);
1342            Double_t pYC = pTC * TMath::Sin(phiC);
1343            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1344            rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1345            rC->SetBgEnergy(0,0);
1346            rC->SetEffArea(jetArea,0);
1347          }
1348        }
1349      }
1350    }// if(fNRandomCones
1351
1352    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1353
1354
1355
1356
1357
1358    if(fAODJetBackgroundOut){
1359      vector<fastjet::PseudoJet> jets2=sortedJets;
1360      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1361      Double_t bkg1=0;
1362      Double_t sigma1=0.;
1363      Double_t meanarea1=0.;
1364      Double_t bkg2=0;
1365      Double_t sigma2=0.;
1366      Double_t meanarea2=0.;
1367
1368      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1369      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1370
1371      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1372      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1373
1374      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1375      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1376      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1377      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1378
1379
1380      //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1381      //exclude 2 leading jets from event
1382      //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R  (Near side to leading track)
1383      //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1384      //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R)  (Away side to leading track)
1385      //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R)   (Orthogonal to leading track)
1386
1387      //Assuming R=0.2 for background clusters
1388
1389      Double_t bkg2Q[4]      = {0.};
1390      Double_t sigma2Q[4]    = {0.};
1391      Double_t meanarea2Q[4] = {0.};
1392
1393      //Get phi of leading track in event
1394      AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1395      Float_t phiLeadingTrack = leading->Phi();
1396
1397      //Quadrant1 - near side
1398      phiMin = phiLeadingTrack - (TMath::Pi()/2./2. - 0.2);
1399      phiMax = phiLeadingTrack + (TMath::Pi()/2./2. - 0.2);
1400      fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1401      clustSeq.get_median_rho_and_sigma(jets2, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1402
1403      //Quadrant2 - orthogonal
1404      phiMin = phiLeadingTrack + TMath::Pi()/2. - (TMath::Pi()/2./2. - 0.2);
1405      phiMax = phiLeadingTrack + TMath::Pi()/2. + (TMath::Pi()/2./2. - 0.2);
1406      fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1407      clustSeq.get_median_rho_and_sigma(jets2, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1408
1409      //Quadrant3 - away side
1410      phiMin = phiLeadingTrack + TMath::Pi() - (TMath::Pi()/2./2. - 0.2);
1411      phiMax = phiLeadingTrack + TMath::Pi() + (TMath::Pi()/2./2. - 0.2);
1412      fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1413      clustSeq.get_median_rho_and_sigma(jets2, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1414
1415      //Quadrant4 - othogonal
1416      phiMin = phiLeadingTrack + TMath::Pi()*3./2. - (TMath::Pi()/2./2. - 0.2);
1417      phiMax = phiLeadingTrack + TMath::Pi()*3./2. + (TMath::Pi()/2./2. - 0.2);
1418      fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1419      clustSeq.get_median_rho_and_sigma(jets2, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1420
1421      fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[0],0.5);
1422      fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[1],1.5);
1423      fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[2],2.5);
1424      fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[3],3.5);
1425
1426      fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[0],0.5);
1427      fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[1],1.5);
1428      fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[2],2.5);
1429      fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[3],3.5);
1430
1431      fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[0],0.5);
1432      fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[1],1.5);
1433      fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[2],2.5);
1434      fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[3],3.5);
1435
1436      fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[0],0.5);
1437      fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[1],1.5);
1438      fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[2],2.5);
1439      fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[3],3.5);
1440
1441    }
1442   }
1443
1444
1445
1446
1447
1448   // fill track information
1449   Int_t nTrackOver = recParticles.GetSize();
1450   // do the same for tracks and jets
1451
1452   if(nTrackOver>0){
1453    TIterator *recIter = recParticles.MakeIterator();
1454    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1455    Float_t pt = tmpRec->Pt();
1456
1457    //    Printf("Leading track p_t %3.3E",pt);
1458    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1459      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1460      while(pt<ptCut&&tmpRec){
1461        nTrackOver--;
1462        tmpRec = (AliVParticle*)(recIter->Next());
1463        if(tmpRec){
1464          pt = tmpRec->Pt();
1465        }
1466      }
1467      if(nTrackOver<=0)break;
1468      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1469    }
1470
1471    recIter->Reset();
1472    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1473    Float_t phi = leading->Phi();
1474    if(phi<0)phi+=TMath::Pi()*2.;
1475    Float_t eta = leading->Eta();
1476    pt  = leading->Pt();
1477    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1478      Float_t tmpPt = tmpRec->Pt();
1479      Float_t tmpEta = tmpRec->Eta();
1480      fh1PtTracksRecIn->Fill(tmpPt);
1481      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1482      if(tmpRec==leading){
1483        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1484        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1485        continue;
1486      }
1487       // correlation
1488      Float_t tmpPhi =  tmpRec->Phi();
1489
1490      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1491      Float_t dPhi = phi - tmpPhi;
1492      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1493      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1494      Float_t dEta = eta - tmpRec->Eta();
1495      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1496      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1497      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1498    }
1499    delete recIter;
1500  }
1501
1502  // find the random jets
1503
1504  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1505
1506  // fill the jet information from random track
1507  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1508  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1509
1510   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1511
1512  // loop over all jets an fill information, first on is the leading jet
1513
1514  Int_t nRecOverRan = inclusiveJetsRan.size();
1515  Int_t nRecRan     = inclusiveJetsRan.size();
1516
1517  if(inclusiveJetsRan.size()>0){
1518    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1519    Float_t pt = leadingJet.Pt();
1520
1521    Int_t iCount = 0;
1522    TLorentzVector vecarearanb;
1523
1524    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1525      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1526       while(pt<ptCut&&iCount<nRecRan){
1527         nRecOverRan--;
1528         iCount++;
1529         if(iCount<nRecRan){
1530           pt = sortedJetsRan[iCount].perp();
1531         }
1532       }
1533       if(nRecOverRan<=0)break;
1534       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1535     }
1536     Float_t phi = leadingJet.Phi();
1537     if(phi<0)phi+=TMath::Pi()*2.;
1538     pt = leadingJet.Pt();
1539
1540     // correlation of leading jet with random tracks
1541
1542     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1543       {
1544         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1545         // correlation
1546         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1547         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1548         Float_t dPhi = phi - tmpPhi;
1549         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1550         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1551         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1552         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1553       }
1554
1555     Int_t nAodOutJetsRan = 0;
1556      AliAODJet *aodOutJetRan = 0;
1557     for(int j = 0; j < nRecRan;j++){
1558       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1559       Float_t tmpPt = tmpRec.Pt();
1560       fh1PtJetsRecInRan->Fill(tmpPt);
1561       // Fill Spectra with constituents
1562       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1563       fh1NConstRecRan->Fill(constituents.size());
1564       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1565       fh2PtNchRan->Fill(nCh,tmpPt);
1566       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1567
1568
1569      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1570        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1571        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1572        aodOutJetRan->GetRefTracks()->Clear();
1573        aodOutJetRan->SetEffArea(arearan,0);
1574        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1575        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1576        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1577
1578      }
1579
1580       // correlation
1581       Float_t tmpPhi =  tmpRec.Phi();
1582       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1583
1584       if(j==0){
1585         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1586         fh1NConstLeadingRecRan->Fill(constituents.size());
1587         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1588         continue;
1589       }
1590     }
1591
1592
1593     if(fAODJetBackgroundOut){
1594      Double_t bkg3=0.;
1595      Double_t sigma3=0.;
1596      Double_t meanarea3=0.;
1597      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1598      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1599      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1600      /*
1601      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1602      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1603      */
1604     }
1605
1606
1607
1608  }
1609
1610
1611  // do the event selection if activated
1612  if(fJetTriggerPtCut>0){
1613    bool select = false;
1614    Float_t minPt = fJetTriggerPtCut;
1615    /*
1616    // hard coded for now ...
1617    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1618    if(cent<10)minPt = 50;
1619    else if(cent<30)minPt = 42;
1620    else if(cent<50)minPt = 28;
1621    else if(cent<80)minPt = 18;
1622    */
1623    float rho = 0;
1624    if(externalBackground)rho = externalBackground->GetBackground(2);
1625    if(fTCAJetsOut){
1626      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1627        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1628        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1629        if(ptSub>=minPt){
1630          select = true;
1631          break;
1632        }
1633      }
1634    }
1635
1636    if(select){
1637      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1638      fh1CentralitySelect->Fill(cent);
1639      fh1ZSelect->Fill(zVtx);
1640      aodH->SetFillAOD(kTRUE);
1641    }
1642  }
1643  if (fDebug > 2){
1644    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1645    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1646    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1647    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1648  }
1649  PostData(1, fHistList);
1650 }
1651
1652 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1653 {
1654   //
1655   // Terminate analysis
1656   //
1657     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1658
1659     if(fMomResH1Fit) delete fMomResH1Fit;
1660     if(fMomResH2Fit) delete fMomResH2Fit;
1661     if(fMomResH3Fit) delete fMomResH3Fit;
1662
1663 }
1664
1665
1666 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1667
1668   //
1669   // get list of tracks/particles for different types
1670   //
1671
1672   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1673
1674   Int_t iCount = 0;
1675   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1676     if(type!=kTrackAODextraonly) {
1677       AliAODEvent *aod = 0;
1678       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1679       else aod = AODEvent();
1680       if(!aod){
1681         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1682         return iCount;
1683       }
1684
1685       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1686         AliAODTrack *tr = aod->GetTrack(it);
1687         Bool_t bGood = false;
1688         if(fFilterType == 0)bGood = true;
1689         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1690         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1691         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1692           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1693           continue;
1694         }
1695         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1696           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1697           continue;
1698         }
1699         if(tr->Pt()<fTrackPtCut){
1700           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1701           continue;
1702         }
1703         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1704         list->Add(tr);
1705         iCount++;
1706       }
1707     }
1708     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1709       AliAODEvent *aod = 0;
1710       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1711       else aod = AODEvent();
1712
1713       if(!aod){
1714         return iCount;
1715       }
1716       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1717       if(!aodExtraTracks)return iCount;
1718       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1719         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1720         if (!track) continue;
1721
1722         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1723         if(!trackAOD)continue;
1724         Bool_t bGood = false;
1725         if(fFilterType == 0)bGood = true;
1726         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1727         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1728         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1729         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1730         if(trackAOD->Pt()<fTrackPtCut) continue;
1731         list->Add(trackAOD);
1732         iCount++;
1733       }
1734     }
1735   }
1736   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1737     AliMCEvent* mcEvent = MCEvent();
1738     if(!mcEvent)return iCount;
1739     // we want to have alivpartilces so use get track
1740     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1741       if(!mcEvent->IsPhysicalPrimary(it))continue;
1742       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1743       if(type == kTrackKineAll){
1744         if(part->Pt()<fTrackPtCut)continue;
1745         list->Add(part);
1746         iCount++;
1747       }
1748       else if(type == kTrackKineCharged){
1749         if(part->Particle()->GetPDG()->Charge()==0)continue;
1750         if(part->Pt()<fTrackPtCut)continue;
1751         list->Add(part);
1752         iCount++;
1753       }
1754     }
1755   }
1756   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1757     AliAODEvent *aod = 0;
1758     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1759     else aod = AODEvent();
1760     if(!aod)return iCount;
1761     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1762     if(!tca)return iCount;
1763     for(int it = 0;it < tca->GetEntriesFast();++it){
1764       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1765       if(!part->IsPhysicalPrimary())continue;
1766       if(type == kTrackAODMCAll){
1767         if(part->Pt()<fTrackPtCut)continue;
1768         list->Add(part);
1769         iCount++;
1770       }
1771       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1772         if(part->Charge()==0)continue;
1773         if(part->Pt()<fTrackPtCut)continue;
1774         if(kTrackAODMCCharged){
1775           list->Add(part);
1776         }
1777         else {
1778           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1779           list->Add(part);
1780         }
1781         iCount++;
1782       }
1783     }
1784   }// AODMCparticle
1785   list->Sort();
1786   return iCount;
1787 }
1788
1789 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1790
1791   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1792
1793   if(!f)return;
1794
1795   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1796   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1797   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1798
1799   SetSmearResolution(kTRUE);
1800   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1801
1802
1803 }
1804
1805 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1806
1807   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1808   if(!f)return;
1809
1810   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1811   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1812   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1813
1814   SetDiceEfficiency(kTRUE);
1815
1816   if(fChangeEfficiencyFraction>0.) {
1817
1818     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1819
1820     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1821       Double_t content = hEffPosGlobSt->GetBinContent(bin);
1822       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1823     }
1824
1825     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1826
1827   }
1828   else
1829     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1830
1831 }
1832
1833 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1834
1835   //
1836   // set mom res profiles
1837   //
1838
1839   if(fMomResH1) delete fMomResH1;
1840   if(fMomResH2) delete fMomResH2;
1841   if(fMomResH3) delete fMomResH3;
1842
1843   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1844   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1845   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1846
1847 }
1848
1849 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1850   //
1851   // set tracking efficiency histos
1852   //
1853
1854   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1855   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1856   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1857 }
1858
1859 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1860
1861   //
1862   // Get smearing on generated momentum
1863   //
1864
1865   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1866
1867   TProfile *fMomRes = 0x0;
1868   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1869   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1870   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1871
1872   if(!fMomRes) {
1873     return 0.;
1874   }
1875
1876
1877   Double_t smear = 0.;
1878
1879   if(pt>20.) {
1880     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1881     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1882     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1883   }
1884   else {
1885
1886     Int_t bin = fMomRes->FindBin(pt);
1887
1888     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1889
1890   }
1891
1892   if(fMomRes) delete fMomRes;
1893
1894   return smear;
1895 }
1896
1897 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1898   //
1899   // Fit linear function on momentum resolution at high pT
1900   //
1901
1902   if(!fMomResH1Fit && fMomResH1) {
1903     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1904     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1905     fMomResH1Fit ->SetRange(5.,100.);
1906   }
1907
1908   if(!fMomResH2Fit && fMomResH2) {
1909     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1910     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1911     fMomResH2Fit ->SetRange(5.,100.);
1912   }
1913
1914   if(!fMomResH3Fit && fMomResH3) {
1915     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1916     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1917     fMomResH3Fit ->SetRange(5.,100.);
1918   }
1919
1920 }
1921
1922 /*
1923 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1924   for(int i = 0; i < particles.GetEntries(); i++){
1925     AliVParticle *vp = (AliVParticle*)particles.At(i);
1926     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1927     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1928     jInp.set_user_index(i);
1929     inputParticles.push_back(jInp);
1930   }
1931
1932   return 0;
1933
1934 }
1935 */