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