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
 343       fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
 344       fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 345       AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
 346
 347       if(fUseBackgroundCalc){
 348         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 349           fAODJetBackgroundOut = new AliAODJetEventBackground();
 350           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 351           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
 352         }
 353       }
 354       // create the branch for the random cones with the same R
 355       TString cName = Form("%sRandomCone",fNonStdBranch.Data());
 356
 357       if(fNRandomCones>0){
 358         if(!AODEvent()->FindListObject(cName.Data())){
 359           fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
 360           fTCARandomConesOut->SetName(cName.Data());
 361           AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
 362         }
 363         // create the branch with the random for the random cones on the random event
 364         cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 365         if(!AODEvent()->FindListObject(cName.Data())){
 366           fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
 367           fTCARandomConesOutRan->SetName(cName.Data());
 368           AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
 369         }
 370       }
 371
 372       if(fNonStdFile.Length()!=0){
 373         //
 374         // case that we have an AOD extension we need to fetch the jets from the extended output
 375         // we identify the extension aod event by looking for the branchname
 376         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 377         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 378         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 379       }
 380     }
 381
 382
 383   if(!fHistList)fHistList = new TList();
 384   fHistList->SetOwner();
 385   PostData(1, fHistList); // post data in any case once
 386
 387   Bool_t oldStatus = TH1::AddDirectoryStatus();
 388   TH1::AddDirectory(kFALSE);
 389
 390   //
 391   //  Histogram
 392
 393   const Int_t nBinPt = 200;
 394   Double_t binLimitsPt[nBinPt+1];
 395   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 396     if(iPt == 0){
 397       binLimitsPt[iPt] = 0.0;
 398     }
 399     else {// 1.0
 400       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 401     }
 402   }
 403
 404   const Int_t nBinPhi = 90;
 405   Double_t binLimitsPhi[nBinPhi+1];
 406   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 407     if(iPhi==0){
 408       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 409     }
 410     else{
 411       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 412     }
 413   }
 414
 415
 416
 417   const Int_t nBinEta = 40;
 418   Double_t binLimitsEta[nBinEta+1];
 419   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 420     if(iEta==0){
 421       binLimitsEta[iEta] = -2.0;
 422     }
 423     else{
 424       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 425     }
 426   }
 427
 428   const int nChMax = 4000;
 429
 430   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 431   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 432
 433   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 434   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 435
 436
 437   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 438   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 439
 440   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 441   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 442   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 443   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 444
 445
 446   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 447   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 448   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 449
 450   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 451   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 452   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 453   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 454   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 455   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 456   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 457   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 458   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 459   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 460
 461   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 462   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
 463   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 464
 465   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 466   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
 467   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 468
 469   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 470   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 471   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 472   //
 473
 474
 475   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 476   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 477   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 478   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 479
 480   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 481   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);
 482   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);
 483   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);
 484
 485
 486
 487   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 488                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 489   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 490                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 491
 492   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 493                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 494   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 495                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 496
 497   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 498                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 499   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 500                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 501
 502
 503
 504   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 505                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 506   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 507                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 508   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 509                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 510   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 511                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 512   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 513                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 514   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 515                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 516
 517   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 518                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 519   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 520                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 521
 522   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 523                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 524   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 525                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 526
 527
 528   if(fNRandomCones>0&&fUseBackgroundCalc){
 529     for(int i = 0;i<3;i++){
 530       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 531       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 532     }
 533   }
 534
 535   for(int i = 0;i < kMaxCent;i++){
 536     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
 537     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
 538     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
 539     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
 540   }
 541
 542   const Int_t saveLevel = 3; // large save level more histos
 543   if(saveLevel>0){
 544     fHistList->Add(fh1Xsec);
 545     fHistList->Add(fh1Trials);
 546
 547     fHistList->Add(fh1NJetsRec);
 548     fHistList->Add(fh1NConstRec);
 549     fHistList->Add(fh1NConstLeadingRec);
 550     fHistList->Add(fh1PtJetsRecIn);
 551     fHistList->Add(fh1PtJetsLeadingRecIn);
 552     fHistList->Add(fh1PtTracksRecIn);
 553     fHistList->Add(fh1PtTracksLeadingRecIn);
 554     fHistList->Add(fh1PtJetConstRec);
 555     fHistList->Add(fh1PtJetConstLeadingRec);
 556     fHistList->Add(fh1NJetsRecRan);
 557     fHistList->Add(fh1NConstRecRan);
 558     fHistList->Add(fh1PtJetsLeadingRecInRan);
 559     fHistList->Add(fh1NConstLeadingRecRan);
 560     fHistList->Add(fh1PtJetsRecInRan);
 561     fHistList->Add(fh1Nch);
 562     fHistList->Add(fh1Centrality);
 563     fHistList->Add(fh1CentralitySelect);
 564     fHistList->Add(fh1CentralityPhySel);
 565     fHistList->Add(fh1Z);
 566     fHistList->Add(fh1ZSelect);
 567     fHistList->Add(fh1ZPhySel);
 568     if(fNRandomCones>0&&fUseBackgroundCalc){
 569       for(int i = 0;i<3;i++){
 570         fHistList->Add(fh1BiARandomCones[i]);
 571         fHistList->Add(fh1BiARandomConesRan[i]);
 572       }
 573     }
 574   for(int i = 0;i < kMaxCent;i++){
 575     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
 576     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
 577     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
 578     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
 579   }
 580
 581     fHistList->Add(fh2NRecJetsPt);
 582     fHistList->Add(fh2NRecTracksPt);
 583     fHistList->Add(fh2NConstPt);
 584     fHistList->Add(fh2NConstLeadingPt);
 585     fHistList->Add(fh2PtNch);
 586     fHistList->Add(fh2PtNchRan);
 587     fHistList->Add(fh2PtNchN);
 588     fHistList->Add(fh2PtNchNRan);
 589     fHistList->Add(fh2JetPhiEta);
 590     fHistList->Add(fh2LeadingJetPhiEta);
 591     fHistList->Add(fh2JetEtaPt);
 592     fHistList->Add(fh2LeadingJetEtaPt);
 593     fHistList->Add(fh2TrackEtaPt);
 594     fHistList->Add(fh2LeadingTrackEtaPt);
 595     fHistList->Add(fh2JetsLeadingPhiEta );
 596     fHistList->Add(fh2JetsLeadingPhiPt);
 597     fHistList->Add(fh2TracksLeadingPhiEta);
 598     fHistList->Add(fh2TracksLeadingPhiPt);
 599     fHistList->Add(fh2TracksLeadingJetPhiPt);
 600     fHistList->Add(fh2JetsLeadingPhiPtW);
 601     fHistList->Add(fh2TracksLeadingPhiPtW);
 602     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 603     fHistList->Add(fh2NRecJetsPtRan);
 604     fHistList->Add(fh2NConstPtRan);
 605     fHistList->Add(fh2NConstLeadingPtRan);
 606     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 607     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 608   }
 609
 610   // =========== Switch on Sumw2 for all histos ===========
 611   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 612     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 613     if (h1){
 614       h1->Sumw2();
 615       continue;
 616     }
 617     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 618     if(hn)hn->Sumw2();
 619   }
 620   TH1::AddDirectory(oldStatus);
 621 }
 622
 623 void AliAnalysisTaskJetCluster::Init()
 624 {
 625   //
 626   // Initialization
 627   //
 628
 629   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 630
 631 }
 632
 633 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 634 {
 635
 636   if(fUseGlobalSelection){
 637     // no selection by the service task, we continue
 638     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 639     PostData(1, fHistList);
 640     return;
 641   }
 642
 643
 644
 645   // handle and reset the output jet branch
 646
 647   if(fTCAJetsOut)fTCAJetsOut->Delete();
 648   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
 649   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
 650   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
 651   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
 652
 653   AliAODJetEventBackground* externalBackground = 0;
 654   if(!externalBackground&&fBackgroundBranch.Length()){
 655     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
 656     if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 657     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 658   }
 659   //
 660   // Execute analysis for current event
 661   //
 662   AliESDEvent *fESD = 0;
 663   if(fUseAODTrackInput){
 664     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 665     if(!fAOD){
 666       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 667       return;
 668     }
 669     // fethc the header
 670   }
 671   else{
 672     //  assume that the AOD is in the general output...
 673     fAOD  = AODEvent();
 674     if(!fAOD){
 675       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 676       return;
 677     }
 678     if(fDebug>0){
 679       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 680     }
 681   }
 682
 683   Bool_t selectEvent =  false;
 684   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 685
 686   Float_t cent = 0;
 687   Float_t zVtx  = 0;
 688   Int_t cenClass = -1;
 689   if(fAOD){
 690     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 691     TString vtxTitle(vtxAOD->GetTitle());
 692     zVtx = vtxAOD->GetZ();
 693
 694     cent = fAOD->GetHeader()->GetCentrality();
 695     if(cent<10)cenClass = 0;
 696     else if(cent<30)cenClass = 1;
 697     else if(cent<50)cenClass = 2;
 698     else if(cent<80)cenClass = 3;
 699     if(physicsSelection){
 700       fh1CentralityPhySel->Fill(cent);
 701       fh1ZPhySel->Fill(zVtx);
 702     }
 703
 704
 705     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 706         Float_t yvtx = vtxAOD->GetY();
 707         Float_t xvtx = vtxAOD->GetX();
 708         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 709         if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
 710           if(physicsSelection){
 711             selectEvent = true;
 712           }
 713         }
 714     }
 715     if(fCentCutUp>0){
 716       if(cent<fCentCutLo||cent>fCentCutUp){
 717         selectEvent = false;
 718       }
 719     }
 720
 721   }
 722   if(!selectEvent){
 723     PostData(1, fHistList);
 724     return;
 725   }
 726   fh1Centrality->Fill(cent);
 727   fh1Z->Fill(zVtx);
 728   fh1Trials->Fill("#sum{ntrials}",1);
 729
 730
 731   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 732
 733   // ==== General variables needed
 734
 735
 736
 737   // we simply fetch the tracks/mc particles as a list of AliVParticles
 738
 739   TList recParticles;
 740   TList genParticles;
 741
 742   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 743   Float_t nCh = recParticles.GetEntries();
 744   fh1Nch->Fill(nCh);
 745   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 746   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 747   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 748
 749   // find the jets....
 750
 751   vector<fastjet::PseudoJet> inputParticlesRec;
 752   vector<fastjet::PseudoJet> inputParticlesRecRan;
 753
 754   // Generate the random cones
 755
 756   AliAODJet vTmpRan(1,0,0,1);
 757   for(int i = 0; i < recParticles.GetEntries(); i++){
 758     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 759     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 760     // we take total momentum here
 761     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 762     jInp.set_user_index(i);
 763     inputParticlesRec.push_back(jInp);
 764
 765     // the randomized input changes eta and phi, but keeps the p_T
 766     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 767       Double_t pT = vp->Pt();
 768       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
 769       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 770
 771       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 772       Double_t pZ = pT/TMath::Tan(theta);
 773
 774       Double_t pX = pT * TMath::Cos(phi);
 775       Double_t pY = pT * TMath::Sin(phi);
 776       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 777       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 778
 779       jInpRan.set_user_index(i);
 780       inputParticlesRecRan.push_back(jInpRan);
 781       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
 782     }
 783
 784     // fill the tref array, only needed when we write out jets
 785     if(fTCAJetsOut){
 786       if(i == 0){
 787         fRef->Delete(); // make sure to delete before placement new...
 788         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
 789       }
 790       fRef->Add(vp);
 791     }
 792   }// recparticles
 793
 794   if(inputParticlesRec.size()==0){
 795     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 796     PostData(1, fHistList);
 797     return;
 798   }
 799
 800   // run fast jet
 801   // employ setters for these...
 802
 803
 804   // now create the object that holds info about ghosts
 805   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
 806     // reduce CPU time...
 807     fGhostArea = 0.5;
 808     fActiveAreaRepeats = 0;
 809   }
 810
 811  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 812   fastjet::AreaType areaType =   fastjet::active_area;
 813   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 814   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 815   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 816   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 817
 818   //range where to compute background
 819   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 820   phiMin = 0;
 821   phiMax = 2*TMath::Pi();
 822   rapMax = fGhostEtamax - fRparam;
 823   rapMin = - fGhostEtamax + fRparam;
 824   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 825
 826
 827
 828   inclusiveJets = clustSeq.inclusive_jets();
 829   sortedJets    = sorted_by_pt(inclusiveJets);
 830
 831   fh1NJetsRec->Fill(sortedJets.size());
 832
 833  // loop over all jets an fill information, first on is the leading jet
 834
 835   Int_t nRecOver = inclusiveJets.size();
 836   Int_t nRec     = inclusiveJets.size();
 837   if(inclusiveJets.size()>0){
 838     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 839     Double_t area = clustSeq.area(sortedJets[0]);
 840     leadingJet.SetEffArea(area,0);
 841     Float_t pt = leadingJet.Pt();
 842     Int_t nAodOutJets = 0;
 843     Int_t nAodOutTracks = 0;
 844     AliAODJet *aodOutJet = 0;
 845
 846     Int_t iCount = 0;
 847     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 848       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 849       while(pt<ptCut&&iCount<nRec){
 850         nRecOver--;
 851         iCount++;
 852         if(iCount<nRec){
 853           pt = sortedJets[iCount].perp();
 854         }
 855       }
 856       if(nRecOver<=0)break;
 857       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 858     }
 859     Float_t phi = leadingJet.Phi();
 860     if(phi<0)phi+=TMath::Pi()*2.;
 861     Float_t eta = leadingJet.Eta();
 862     Float_t pTback = 0;
 863     if(externalBackground){
 864       // carefull has to be filled in a task before
 865       // todo, ReArrange to the botom
 866       pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
 867     }
 868     pt = leadingJet.Pt() - pTback;
 869     // correlation of leading jet with tracks
 870     TIterator *recIter = recParticles.MakeIterator();
 871     recIter->Reset();
 872     AliVParticle *tmpRecTrack = 0;
 873     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 874       Float_t tmpPt = tmpRecTrack->Pt();
 875       // correlation
 876       Float_t tmpPhi =  tmpRecTrack->Phi();
 877       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 878       Float_t dPhi = phi - tmpPhi;
 879       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 880       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 881       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 882       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 883       if(cenClass>=0){
 884         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
 885         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
 886       }
 887
 888     }
 889
 890
 891     TLorentzVector vecareab;
 892     for(int j = 0; j < nRec;j++){
 893       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 894       aodOutJet = 0;
 895       nAodOutTracks = 0;
 896       Float_t tmpPt = tmpRec.Pt();
 897
 898       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
 899         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
 900         aodOutJet->GetRefTracks()->Clear();
 901         Double_t area1 = clustSeq.area(sortedJets[j]);
 902         aodOutJet->SetEffArea(area1,0);
 903         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
 904         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
 905         aodOutJet->SetVectorAreaCharged(&vecareab);
 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 = fTrackEtaWindow - 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        aodOutJetRan->GetRefTracks()->Clear();
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       // correlation
1262       Float_t tmpPhi =  tmpRec.Phi();
1263       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1264
1265       if(j==0){
1266         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1267         fh1NConstLeadingRecRan->Fill(constituents.size());
1268         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1269         continue;
1270       }
1271     }
1272
1273
1274     if(fAODJetBackgroundOut){
1275      Double_t bkg3=0.;
1276      Double_t sigma3=0.;
1277      Double_t meanarea3=0.;
1278      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1279      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1280      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1281      /*
1282      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1283      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1284      */
1285     }
1286
1287
1288
1289  }
1290
1291
1292  // do the event selection if activated
1293  if(fJetTriggerPtCut>0){
1294    bool select = false;
1295    Float_t minPt = fJetTriggerPtCut;
1296    /*
1297    // hard coded for now ...
1298    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1299    if(cent<10)minPt = 50;
1300    else if(cent<30)minPt = 42;
1301    else if(cent<50)minPt = 28;
1302    else if(cent<80)minPt = 18;
1303    */
1304    float rho = 0;
1305    if(externalBackground)rho = externalBackground->GetBackground(2);
1306    if(fTCAJetsOut){
1307      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1308        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1309        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1310        if(ptSub>=minPt){
1311          select = true;
1312          break;
1313        }
1314      }
1315    }
1316
1317    if(select){
1318      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1319      fh1CentralitySelect->Fill(cent);
1320      fh1ZSelect->Fill(zVtx);
1321      aodH->SetFillAOD(kTRUE);
1322    }
1323  }
1324
1325  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1326  PostData(1, fHistList);
1327 }
1328
1329 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1330 {
1331 // Terminate analysis
1332 //
1333     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1334 }
1335
1336
1337 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1338
1339   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1340
1341   Int_t iCount = 0;
1342   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1343     if(type!=kTrackAODextraonly) {
1344       AliAODEvent *aod = 0;
1345       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1346       else aod = AODEvent();
1347       if(!aod){
1348         return iCount;
1349       }
1350       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1351         AliAODTrack *tr = aod->GetTrack(it);
1352         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1353         if(TMath::Abs(tr->Eta())>fTrackEtaWindow)continue;
1354         if(tr->Pt()<fTrackPtCut)continue;
1355         list->Add(tr);
1356         iCount++;
1357       }
1358     }
1359     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1360       AliAODEvent *aod = 0;
1361       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1362       else aod = AODEvent();
1363
1364       if(!aod){
1365         return iCount;
1366       }
1367       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1368       if(!aodExtraTracks)return iCount;
1369       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1370         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1371         if (!track) continue;
1372
1373         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1374         if(!trackAOD)continue;
1375         if(trackAOD->Pt()<fTrackPtCut) continue;
1376         list->Add(trackAOD);
1377         iCount++;
1378       }
1379     }
1380   }
1381   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1382     AliMCEvent* mcEvent = MCEvent();
1383     if(!mcEvent)return iCount;
1384     // we want to have alivpartilces so use get track
1385     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1386       if(!mcEvent->IsPhysicalPrimary(it))continue;
1387       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1388       if(type == kTrackKineAll){
1389         if(part->Pt()<fTrackPtCut)continue;
1390         list->Add(part);
1391         iCount++;
1392       }
1393       else if(type == kTrackKineCharged){
1394         if(part->Particle()->GetPDG()->Charge()==0)continue;
1395         if(part->Pt()<fTrackPtCut)continue;
1396         list->Add(part);
1397         iCount++;
1398       }
1399     }
1400   }
1401   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1402     AliAODEvent *aod = 0;
1403     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1404     else aod = AODEvent();
1405     if(!aod)return iCount;
1406     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1407     if(!tca)return iCount;
1408     for(int it = 0;it < tca->GetEntriesFast();++it){
1409       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1410       if(!part->IsPhysicalPrimary())continue;
1411       if(type == kTrackAODMCAll){
1412         if(part->Pt()<fTrackPtCut)continue;
1413         list->Add(part);
1414         iCount++;
1415       }
1416       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1417         if(part->Charge()==0)continue;
1418         if(part->Pt()<fTrackPtCut)continue;
1419         if(kTrackAODMCCharged){
1420           list->Add(part);
1421         }
1422         else {
1423           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1424           list->Add(part);
1425         }
1426         iCount++;
1427       }
1428     }
1429   }// AODMCparticle
1430   list->Sort();
1431   return iCount;
1432 }
1433
1434 /*
1435 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1436   for(int i = 0; i < particles.GetEntries(); i++){
1437     AliVParticle *vp = (AliVParticle*)particles.At(i);
1438     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1439     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1440     jInp.set_user_index(i);
1441     inputParticles.push_back(jInp);
1442   }
1443
1444   return 0;
1445
1446 }
1447 */