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