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