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