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