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