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
1142       AliAODTrack *aodT = 0;
1143       if(partLead){
1144         aodT = dynamic_cast<AliAODTrack*>(partLead);
1145         if(aodT){
1146           if(aodT->TestFilterBit(fFilterMaskBestPt)){
1147             aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1148           }
1149         }
1150       }
1151
1152      // correlation
1153      Float_t tmpPhi =  tmpRec.Phi();
1154      Float_t tmpEta =  tmpRec.Eta();
1155      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1156      if(j==0){
1157        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1158        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1159        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1160        fh1NConstLeadingRec->Fill(constituents.size());
1161        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1162        continue;
1163      }
1164      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1165      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1166      Float_t dPhi = phi - tmpPhi;
1167      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1168      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1169      Float_t dEta = eta - tmpRec.Eta();
1170      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1171      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1172      if(cenClass>=0){
1173        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1174        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1175      }
1176      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1177     }// loop over reconstructed jets
1178    delete recIter;
1179
1180
1181
1182    // Add the random cones...
1183    if(fNRandomCones>0&&fTCARandomConesOut){
1184      // create a random jet within the acceptance
1185      Double_t etaMax = fTrackEtaWindow - fRparam;
1186      Int_t nCone = 0;
1187      Int_t nConeRan = 0;
1188      Double_t pTC = 1; // small number
1189      for(int ir = 0;ir < fNRandomCones;ir++){
1190        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1191        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1192        // massless jet
1193        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1194        Double_t pZC = pTC/TMath::Tan(thetaC);
1195        Double_t pXC = pTC * TMath::Cos(phiC);
1196        Double_t pYC = pTC * TMath::Sin(phiC);
1197        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1198        AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1199        bool skip = false;
1200        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1201          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1202          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1203            skip = true;
1204            break;
1205          }
1206        }
1207        // test for overlap with previous cones to avoid double counting
1208        for(int iic = 0;iic<ir;iic++){
1209          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1210          if(iicone){
1211            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1212              skip = true;
1213              break;
1214            }
1215          }
1216        }
1217        if(skip)continue;
1218        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1219        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1220        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1221      }// loop over random cones creation
1222
1223
1224      // loop over the reconstructed particles and add up the pT in the random cones
1225      // maybe better to loop over randomized particles not in the real jets...
1226      // but this by definition brings dow average energy in the whole  event
1227      AliAODJet vTmpRanR(1,0,0,1);
1228      for(int i = 0; i < recParticles.GetEntries(); i++){
1229        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1230        if(fTCARandomConesOut){
1231          for(int ir = 0;ir < fNRandomCones;ir++){
1232            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1233            if(jC&&jC->DeltaR(vp)<fRparam){
1234              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1235              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1236            }
1237          }
1238        }// add up energy in cone
1239
1240        // the randomized input changes eta and phi, but keeps the p_T
1241        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1242          Double_t pTR = vp->Pt();
1243          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1244          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1245
1246          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1247          Double_t pZR = pTR/TMath::Tan(thetaR);
1248
1249          Double_t pXR = pTR * TMath::Cos(phiR);
1250          Double_t pYR = pTR * TMath::Sin(phiR);
1251          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
1252          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1253          if(fTCARandomConesOutRan){
1254            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1255              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1256              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1257                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1258                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1259              }
1260            }
1261          }
1262        }
1263      }// loop over recparticles
1264
1265      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1266      if(fTCARandomConesOut){
1267        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1268          // rescale the momntum vectors for the random cones
1269
1270          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1271          if(rC){
1272            Double_t etaC = rC->Eta();
1273            Double_t phiC = rC->Phi();
1274            // massless jet, unit vector
1275            pTC = rC->ChargedBgEnergy();
1276            if(pTC<=0)pTC = 0.001; // for almost empty events
1277            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1278            Double_t pZC = pTC/TMath::Tan(thetaC);
1279            Double_t pXC = pTC * TMath::Cos(phiC);
1280            Double_t pYC = pTC * TMath::Sin(phiC);
1281            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1282            rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1283            rC->SetBgEnergy(0,0);
1284            rC->SetEffArea(jetArea,0);
1285          }
1286        }
1287      }
1288      if(fTCARandomConesOutRan){
1289        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1290          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1291          // same wit random
1292          if(rC){
1293            Double_t etaC = rC->Eta();
1294            Double_t phiC = rC->Phi();
1295            // massless jet, unit vector
1296            pTC = rC->ChargedBgEnergy();
1297            if(pTC<=0)pTC = 0.001;// for almost empty events
1298            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1299            Double_t pZC = pTC/TMath::Tan(thetaC);
1300            Double_t pXC = pTC * TMath::Cos(phiC);
1301            Double_t pYC = pTC * TMath::Sin(phiC);
1302            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1303            rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1304            rC->SetBgEnergy(0,0);
1305            rC->SetEffArea(jetArea,0);
1306          }
1307        }
1308      }
1309    }// if(fNRandomCones
1310
1311    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1312
1313
1314
1315
1316
1317    if(fAODJetBackgroundOut){
1318      vector<fastjet::PseudoJet> jets2=sortedJets;
1319      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1320      Double_t bkg1=0;
1321      Double_t sigma1=0.;
1322      Double_t meanarea1=0.;
1323      Double_t bkg2=0;
1324      Double_t sigma2=0.;
1325      Double_t meanarea2=0.;
1326
1327      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1328      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1329
1330      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1331      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1332
1333      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1334      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1335      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1336      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1337
1338    }
1339   }
1340
1341
1342
1343
1344
1345   // fill track information
1346   Int_t nTrackOver = recParticles.GetSize();
1347   // do the same for tracks and jets
1348
1349   if(nTrackOver>0){
1350    TIterator *recIter = recParticles.MakeIterator();
1351    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1352    Float_t pt = tmpRec->Pt();
1353
1354    //    Printf("Leading track p_t %3.3E",pt);
1355    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1356      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1357      while(pt<ptCut&&tmpRec){
1358        nTrackOver--;
1359        tmpRec = (AliVParticle*)(recIter->Next());
1360        if(tmpRec){
1361          pt = tmpRec->Pt();
1362        }
1363      }
1364      if(nTrackOver<=0)break;
1365      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1366    }
1367
1368    recIter->Reset();
1369    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1370    Float_t phi = leading->Phi();
1371    if(phi<0)phi+=TMath::Pi()*2.;
1372    Float_t eta = leading->Eta();
1373    pt  = leading->Pt();
1374    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1375      Float_t tmpPt = tmpRec->Pt();
1376      Float_t tmpEta = tmpRec->Eta();
1377      fh1PtTracksRecIn->Fill(tmpPt);
1378      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1379      if(tmpRec==leading){
1380        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1381        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1382        continue;
1383      }
1384       // correlation
1385      Float_t tmpPhi =  tmpRec->Phi();
1386
1387      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1388      Float_t dPhi = phi - tmpPhi;
1389      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1390      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1391      Float_t dEta = eta - tmpRec->Eta();
1392      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1393      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1394      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1395    }
1396    delete recIter;
1397  }
1398
1399  // find the random jets
1400
1401  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1402
1403  // fill the jet information from random track
1404  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1405  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1406
1407   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1408
1409  // loop over all jets an fill information, first on is the leading jet
1410
1411  Int_t nRecOverRan = inclusiveJetsRan.size();
1412  Int_t nRecRan     = inclusiveJetsRan.size();
1413
1414  if(inclusiveJetsRan.size()>0){
1415    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1416    Float_t pt = leadingJet.Pt();
1417
1418    Int_t iCount = 0;
1419    TLorentzVector vecarearanb;
1420
1421    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1422      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1423       while(pt<ptCut&&iCount<nRecRan){
1424         nRecOverRan--;
1425         iCount++;
1426         if(iCount<nRecRan){
1427           pt = sortedJetsRan[iCount].perp();
1428         }
1429       }
1430       if(nRecOverRan<=0)break;
1431       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1432     }
1433     Float_t phi = leadingJet.Phi();
1434     if(phi<0)phi+=TMath::Pi()*2.;
1435     pt = leadingJet.Pt();
1436
1437     // correlation of leading jet with random tracks
1438
1439     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1440       {
1441         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1442         // correlation
1443         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1444         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1445         Float_t dPhi = phi - tmpPhi;
1446         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1447         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1448         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1449         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1450       }
1451
1452     Int_t nAodOutJetsRan = 0;
1453      AliAODJet *aodOutJetRan = 0;
1454     for(int j = 0; j < nRecRan;j++){
1455       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1456       Float_t tmpPt = tmpRec.Pt();
1457       fh1PtJetsRecInRan->Fill(tmpPt);
1458       // Fill Spectra with constituents
1459       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1460       fh1NConstRecRan->Fill(constituents.size());
1461       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1462       fh2PtNchRan->Fill(nCh,tmpPt);
1463       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1464
1465
1466      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1467        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1468        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1469        aodOutJetRan->GetRefTracks()->Clear();
1470        aodOutJetRan->SetEffArea(arearan,0);
1471        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1472        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1473        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1474
1475      }
1476
1477       // correlation
1478       Float_t tmpPhi =  tmpRec.Phi();
1479       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1480
1481       if(j==0){
1482         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1483         fh1NConstLeadingRecRan->Fill(constituents.size());
1484         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1485         continue;
1486       }
1487     }
1488
1489
1490     if(fAODJetBackgroundOut){
1491      Double_t bkg3=0.;
1492      Double_t sigma3=0.;
1493      Double_t meanarea3=0.;
1494      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1495      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1496      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1497      /*
1498      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1499      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1500      */
1501     }
1502
1503
1504
1505  }
1506
1507
1508  // do the event selection if activated
1509  if(fJetTriggerPtCut>0){
1510    bool select = false;
1511    Float_t minPt = fJetTriggerPtCut;
1512    /*
1513    // hard coded for now ...
1514    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1515    if(cent<10)minPt = 50;
1516    else if(cent<30)minPt = 42;
1517    else if(cent<50)minPt = 28;
1518    else if(cent<80)minPt = 18;
1519    */
1520    float rho = 0;
1521    if(externalBackground)rho = externalBackground->GetBackground(2);
1522    if(fTCAJetsOut){
1523      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1524        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1525        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1526        if(ptSub>=minPt){
1527          select = true;
1528          break;
1529        }
1530      }
1531    }
1532
1533    if(select){
1534      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1535      fh1CentralitySelect->Fill(cent);
1536      fh1ZSelect->Fill(zVtx);
1537      aodH->SetFillAOD(kTRUE);
1538    }
1539  }
1540  if (fDebug > 2){
1541    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1542    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1543    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1544    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1545  }
1546  PostData(1, fHistList);
1547 }
1548
1549 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1550 {
1551   //
1552   // Terminate analysis
1553   //
1554     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1555
1556     if(fMomResH1Fit) delete fMomResH1Fit;
1557     if(fMomResH2Fit) delete fMomResH2Fit;
1558     if(fMomResH3Fit) delete fMomResH3Fit;
1559
1560 }
1561
1562
1563 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1564
1565   //
1566   // get list of tracks/particles for different types
1567   //
1568
1569   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1570
1571   Int_t iCount = 0;
1572   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1573     if(type!=kTrackAODextraonly) {
1574       AliAODEvent *aod = 0;
1575       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1576       else aod = AODEvent();
1577       if(!aod){
1578         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1579         return iCount;
1580       }
1581
1582       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1583         AliAODTrack *tr = aod->GetTrack(it);
1584         Bool_t bGood = false;
1585         if(fFilterType == 0)bGood = true;
1586         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1587         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1588         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1589           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1590           continue;
1591         }
1592         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1593           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1594           continue;
1595         }
1596         if(tr->Pt()<fTrackPtCut){
1597           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1598           continue;
1599         }
1600         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1601         list->Add(tr);
1602         iCount++;
1603       }
1604     }
1605     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1606       AliAODEvent *aod = 0;
1607       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1608       else aod = AODEvent();
1609
1610       if(!aod){
1611         return iCount;
1612       }
1613       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1614       if(!aodExtraTracks)return iCount;
1615       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1616         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1617         if (!track) continue;
1618
1619         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1620         if(!trackAOD)continue;
1621         Bool_t bGood = false;
1622         if(fFilterType == 0)bGood = true;
1623         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1624         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1625         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1626         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1627         if(trackAOD->Pt()<fTrackPtCut) continue;
1628         list->Add(trackAOD);
1629         iCount++;
1630       }
1631     }
1632   }
1633   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1634     AliMCEvent* mcEvent = MCEvent();
1635     if(!mcEvent)return iCount;
1636     // we want to have alivpartilces so use get track
1637     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1638       if(!mcEvent->IsPhysicalPrimary(it))continue;
1639       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1640       if(type == kTrackKineAll){
1641         if(part->Pt()<fTrackPtCut)continue;
1642         list->Add(part);
1643         iCount++;
1644       }
1645       else if(type == kTrackKineCharged){
1646         if(part->Particle()->GetPDG()->Charge()==0)continue;
1647         if(part->Pt()<fTrackPtCut)continue;
1648         list->Add(part);
1649         iCount++;
1650       }
1651     }
1652   }
1653   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1654     AliAODEvent *aod = 0;
1655     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1656     else aod = AODEvent();
1657     if(!aod)return iCount;
1658     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1659     if(!tca)return iCount;
1660     for(int it = 0;it < tca->GetEntriesFast();++it){
1661       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1662       if(!part->IsPhysicalPrimary())continue;
1663       if(type == kTrackAODMCAll){
1664         if(part->Pt()<fTrackPtCut)continue;
1665         list->Add(part);
1666         iCount++;
1667       }
1668       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1669         if(part->Charge()==0)continue;
1670         if(part->Pt()<fTrackPtCut)continue;
1671         if(kTrackAODMCCharged){
1672           list->Add(part);
1673         }
1674         else {
1675           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1676           list->Add(part);
1677         }
1678         iCount++;
1679       }
1680     }
1681   }// AODMCparticle
1682   list->Sort();
1683   return iCount;
1684 }
1685
1686 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1687
1688   //
1689   // set mom res profiles
1690   //
1691
1692   fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
1693   fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
1694   fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
1695 }
1696
1697 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1698   //
1699   // set tracking efficiency histos
1700   //
1701
1702   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1703   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1704   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1705 }
1706
1707 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1708
1709   //
1710   // Get smearing on generated momentum
1711   //
1712
1713   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1714
1715   TProfile *fMomRes = 0x0;
1716   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1717   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1718   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1719
1720   if(!fMomRes) {
1721     return 0.;
1722   }
1723
1724
1725   Double_t smear = 0.;
1726
1727   if(pt>20.) {
1728     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1729     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1730     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1731   }
1732   else {
1733
1734     Int_t bin = fMomRes->FindBin(pt);
1735
1736     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1737
1738   }
1739
1740   if(fMomRes) delete fMomRes;
1741
1742   return smear;
1743 }
1744
1745 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1746   //
1747   // Fit linear function on momentum resolution at high pT
1748   //
1749
1750   if(!fMomResH1Fit && fMomResH1) {
1751     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1752     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1753     fMomResH1Fit ->SetRange(5.,100.);
1754   }
1755
1756   if(!fMomResH2Fit && fMomResH2) {
1757     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1758     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1759     fMomResH2Fit ->SetRange(5.,100.);
1760   }
1761
1762   if(!fMomResH3Fit && fMomResH3) {
1763     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1764     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1765     fMomResH3Fit ->SetRange(5.,100.);
1766   }
1767
1768 }
1769
1770 /*
1771 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1772   for(int i = 0; i < particles.GetEntries(); i++){
1773     AliVParticle *vp = (AliVParticle*)particles.At(i);
1774     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1775     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1776     jInp.set_user_index(i);
1777     inputParticles.push_back(jInp);
1778   }
1779
1780   return 0;
1781
1782 }
1783 */