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 identifay 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(2)*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
 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       }
 917
 918
 919       Float_t tmpPtBack = 0;
 920       if(externalBackground){
 921         // carefull has to be filled in a task before
 922        // todo, ReArrange to the botom
 923         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
 924       }
 925       tmpPt = tmpPt - tmpPtBack;
 926       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
 927
 928       fh1PtJetsRecIn->Fill(tmpPt);
 929       // Fill Spectra with constituents
 930       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 931
 932       fh1NConstRec->Fill(constituents.size());
 933       fh2PtNch->Fill(nCh,tmpPt);
 934       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 935       fh2NConstPt->Fill(tmpPt,constituents.size());
 936       // loop over constiutents and fill spectrum
 937
 938       // Add the jet information and the track references to the output AOD
 939
 940
 941       if(fNRandomCones>0){
 942         // create a random jet within the acceptance
 943         Double_t etaMax = 0.8 - fRparam;
 944         Int_t nCone = 0;
 945         Int_t nConeRan = 0;
 946         Double_t pTC = 1; // small number
 947         for(int ir = 0;ir < fNRandomCones;ir++){
 948           Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 949           Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 950           // massless jet
 951           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 952           Double_t pZC = pTC/TMath::Tan(thetaC);
 953           Double_t pXC = pTC * TMath::Cos(phiC);
 954           Double_t pYC = pTC * TMath::Sin(phiC);
 955           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 956           AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 957           bool skip = false;
 958           for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
 959             AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
 960             if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
 961               skip = true;
 962               break;
 963             }
 964           }
 965           // test for overlap with previous cones to avoid double counting
 966           for(int iic = 0;iic<ir;iic++){
 967             AliAODJet *iicone = (AliAODJet*)rConeArray->At(iic);
 968             if(iicone){
 969               if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
 970                 skip = true;
 971                 break;
 972               }
 973             }
 974           }
 975           if(skip)continue;
 976           tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
 977           if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
 978           if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
 979         }// random cones
 980
 981
 982         // loop over the reconstructed particles and add up the pT in the random cones
 983         // maybe better to loop over randomized particles not in the real jets...
 984         // but this by definition brings dow average energy in the whole  event
 985         AliAODJet vTmpRanR(1,0,0,1);
 986         for(int i = 0; i < recParticles.GetEntries(); i++){
 987           AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 988           if(rConeArray){
 989             for(int ir = 0;ir < fNRandomCones;ir++){
 990               AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
 991               if(jC&&jC->DeltaR(vp)<fRparam){
 992                 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
 993               }
 994             }
 995           }// add up energy in cone
 996
 997           // the randomized input changes eta and phi, but keeps the p_T
 998           if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 999             Double_t pTR = vp->Pt();
1000             Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
1001             Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1002
1003             Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1004             Double_t pZR = pTR/TMath::Tan(thetaR);
1005
1006             Double_t pXR = pTR * TMath::Cos(phiR);
1007             Double_t pYR = pTR * TMath::Sin(phiR);
1008             Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
1009             vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1010             if(rConeArrayRan){
1011               for(int ir = 0;ir < fNRandomCones;ir++){
1012                 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
1013                 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1014                   jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1015                 }
1016               }
1017             }
1018           }
1019         }// loop over recparticles
1020
1021         Float_t jetArea = fRparam*fRparam*TMath::Pi();
1022         for(int ir = 0;ir < fNRandomCones;ir++){
1023           // rescale the momntum vectors for the random cones
1024           if(!rConeArray)continue;
1025           AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1026           if(rC){
1027             Double_t etaC = rC->Eta();
1028             Double_t phiC = rC->Phi();
1029             // massless jet, unit vector
1030             pTC = rC->ChargedBgEnergy();
1031             if(pTC<=0)pTC = 0.1; // for almost empty events
1032             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1033             Double_t pZC = pTC/TMath::Tan(thetaC);
1034             Double_t pXC = pTC * TMath::Cos(phiC);
1035             Double_t pYC = pTC * TMath::Sin(phiC);
1036             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1037             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1038             rC->SetBgEnergy(0,0);
1039             rC->SetEffArea(jetArea,0);
1040           }
1041           rC = (AliAODJet*)rConeArrayRan->At(ir);
1042           // same wit random
1043           if(rC){
1044             Double_t etaC = rC->Eta();
1045             Double_t phiC = rC->Phi();
1046             // massless jet, unit vector
1047             pTC = rC->ChargedBgEnergy();
1048             if(pTC<=0)pTC = 0.1;// for almost empty events
1049             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1050             Double_t pZC = pTC/TMath::Tan(thetaC);
1051             Double_t pXC = pTC * TMath::Cos(phiC);
1052             Double_t pYC = pTC * TMath::Sin(phiC);
1053             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1054             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1055             rC->SetBgEnergy(0,0);
1056             rC->SetEffArea(jetArea,0);
1057           }
1058         }
1059       }// if(fUseBackgroundCalc
1060
1061      for(unsigned int ic = 0; ic < constituents.size();ic++){
1062        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1063        fh1PtJetConstRec->Fill(part->Pt());
1064        if(aodOutJet){
1065          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1066        }
1067        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1068      }
1069
1070      // correlation
1071      Float_t tmpPhi =  tmpRec.Phi();
1072      Float_t tmpEta =  tmpRec.Eta();
1073      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1074      if(j==0){
1075        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1076        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1077        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1078        fh1NConstLeadingRec->Fill(constituents.size());
1079        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1080        continue;
1081      }
1082      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1083      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1084      Float_t dPhi = phi - tmpPhi;
1085      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1086      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1087      Float_t dEta = eta - tmpRec.Eta();
1088      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1089      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1090      if(cenClass>=0){
1091        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1092        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1093      }
1094      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1095     }
1096    delete recIter;
1097
1098    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1099
1100
1101    if(evBkg){
1102      vector<fastjet::PseudoJet> jets2=sortedJets;
1103      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1104      Double_t bkg1=0;
1105      Double_t sigma1=0.;
1106      Double_t meanarea1=0.;
1107      Double_t bkg2=0;
1108      Double_t sigma2=0.;
1109      Double_t meanarea2=0.;
1110
1111      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
1112      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1113
1114      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1115      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1116
1117      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1118      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1119      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1120      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1121
1122    }
1123   }
1124
1125
1126
1127
1128
1129   // fill track information
1130   Int_t nTrackOver = recParticles.GetSize();
1131   // do the same for tracks and jets
1132
1133   if(nTrackOver>0){
1134    TIterator *recIter = recParticles.MakeIterator();
1135    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1136    Float_t pt = tmpRec->Pt();
1137
1138    //    Printf("Leading track p_t %3.3E",pt);
1139    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1140      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1141      while(pt<ptCut&&tmpRec){
1142        nTrackOver--;
1143        tmpRec = (AliVParticle*)(recIter->Next());
1144        if(tmpRec){
1145          pt = tmpRec->Pt();
1146        }
1147      }
1148      if(nTrackOver<=0)break;
1149      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1150    }
1151
1152    recIter->Reset();
1153    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1154    Float_t phi = leading->Phi();
1155    if(phi<0)phi+=TMath::Pi()*2.;
1156    Float_t eta = leading->Eta();
1157    pt  = leading->Pt();
1158    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1159      Float_t tmpPt = tmpRec->Pt();
1160      Float_t tmpEta = tmpRec->Eta();
1161      fh1PtTracksRecIn->Fill(tmpPt);
1162      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1163      if(tmpRec==leading){
1164        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1165        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1166        continue;
1167      }
1168       // correlation
1169      Float_t tmpPhi =  tmpRec->Phi();
1170
1171      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1172      Float_t dPhi = phi - tmpPhi;
1173      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1174      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1175      Float_t dEta = eta - tmpRec->Eta();
1176      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1177      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1178      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1179    }
1180    delete recIter;
1181  }
1182
1183  // find the random jets
1184  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1185  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1186
1187  inclusiveJetsRan = clustSeqRan.inclusive_jets();
1188  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1189
1190  // fill the jet information from random track
1191
1192   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1193
1194  // loop over all jets an fill information, first on is the leading jet
1195
1196  Int_t nRecOverRan = inclusiveJetsRan.size();
1197  Int_t nRecRan     = inclusiveJetsRan.size();
1198  if(inclusiveJetsRan.size()>0){
1199    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1200    Float_t pt = leadingJet.Pt();
1201
1202    Int_t iCount = 0;
1203    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1204      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1205       while(pt<ptCut&&iCount<nRecRan){
1206         nRecOverRan--;
1207         iCount++;
1208         if(iCount<nRecRan){
1209           pt = sortedJetsRan[iCount].perp();
1210         }
1211       }
1212       if(nRecOverRan<=0)break;
1213       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1214     }
1215     Float_t phi = leadingJet.Phi();
1216     if(phi<0)phi+=TMath::Pi()*2.;
1217     pt = leadingJet.Pt();
1218
1219     // correlation of leading jet with random tracks
1220
1221     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1222       {
1223         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1224         // correlation
1225         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1226         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1227         Float_t dPhi = phi - tmpPhi;
1228         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1229         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1230         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1231         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1232       }
1233
1234     Int_t nAodOutJetsRan = 0;
1235      AliAODJet *aodOutJetRan = 0;
1236     for(int j = 0; j < nRecRan;j++){
1237       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1238       Float_t tmpPt = tmpRec.Pt();
1239       fh1PtJetsRecInRan->Fill(tmpPt);
1240       // Fill Spectra with constituents
1241       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1242       fh1NConstRecRan->Fill(constituents.size());
1243       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1244       fh2PtNchRan->Fill(nCh,tmpPt);
1245       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1246
1247
1248      if(tmpPt>fJetOutputMinPt&&jarrayran){
1249        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1250        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1251
1252        aodOutJetRan->SetEffArea(arearan,0);    }
1253
1254
1255
1256
1257      for(unsigned int ic = 0; ic < constituents.size();ic++){
1258        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1259        // fh1PtJetConstRec->Fill(part->Pt());
1260        if(aodOutJetRan){
1261          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1262        }
1263      }
1264
1265
1266
1267
1268       // correlation
1269       Float_t tmpPhi =  tmpRec.Phi();
1270       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1271
1272       if(j==0){
1273         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1274         fh1NConstLeadingRecRan->Fill(constituents.size());
1275         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1276         continue;
1277       }
1278     }
1279
1280
1281     if(evBkg){
1282      Double_t bkg3=0.;
1283      Double_t sigma3=0.;
1284      Double_t meanarea3=0.;
1285      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1286      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1287      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1288      /*
1289      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1290      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1291      */
1292     }
1293
1294
1295
1296  }
1297
1298
1299  // do the event selection if activated
1300  if(fJetTriggerPtCut>0){
1301    bool select = false;
1302    Float_t minPt = fJetTriggerPtCut;
1303    /*
1304    // hard coded for now ...
1305    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1306    if(cent<10)minPt = 50;
1307    else if(cent<30)minPt = 42;
1308    else if(cent<50)minPt = 28;
1309    else if(cent<80)minPt = 18;
1310    */
1311    float rho = 0;
1312    if(externalBackground)rho = externalBackground->GetBackground(2);
1313    if(jarray){
1314      for(int i = 0;i < jarray->GetEntriesFast();i++){
1315        AliAODJet *jet = (AliAODJet*)jarray->At(i);
1316        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1317        if(ptSub>=minPt){
1318          select = true;
1319          break;
1320        }
1321      }
1322    }
1323
1324    if(select){
1325      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1326      fh1CentralitySelect->Fill(cent);
1327      fh1ZSelect->Fill(zVtx);
1328      aodH->SetFillAOD(kTRUE);
1329    }
1330  }
1331
1332  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1333  PostData(1, fHistList);
1334 }
1335
1336 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1337 {
1338 // Terminate analysis
1339 //
1340     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1341 }
1342
1343
1344 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1345
1346   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1347
1348   Int_t iCount = 0;
1349   if(type==kTrackAOD){
1350     AliAODEvent *aod = 0;
1351     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1352     else aod = AODEvent();
1353     if(!aod){
1354       return iCount;
1355     }
1356     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1357       AliAODTrack *tr = aod->GetTrack(it);
1358       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1359       if(TMath::Abs(tr->Eta())>0.9)continue;
1360       if(tr->Pt()<fTrackPtCut)continue;
1361       list->Add(tr);
1362       iCount++;
1363     }
1364   }
1365   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1366     AliMCEvent* mcEvent = MCEvent();
1367     if(!mcEvent)return iCount;
1368     // we want to have alivpartilces so use get track
1369     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1370       if(!mcEvent->IsPhysicalPrimary(it))continue;
1371       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1372       if(type == kTrackKineAll){
1373         if(part->Pt()<fTrackPtCut)continue;
1374         list->Add(part);
1375         iCount++;
1376       }
1377       else if(type == kTrackKineCharged){
1378         if(part->Particle()->GetPDG()->Charge()==0)continue;
1379         if(part->Pt()<fTrackPtCut)continue;
1380         list->Add(part);
1381         iCount++;
1382       }
1383     }
1384   }
1385   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1386     AliAODEvent *aod = 0;
1387     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1388     else aod = AODEvent();
1389     if(!aod)return iCount;
1390     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1391     if(!tca)return iCount;
1392     for(int it = 0;it < tca->GetEntriesFast();++it){
1393       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1394       if(!part->IsPhysicalPrimary())continue;
1395       if(type == kTrackAODMCAll){
1396         if(part->Pt()<fTrackPtCut)continue;
1397         list->Add(part);
1398         iCount++;
1399       }
1400       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1401         if(part->Charge()==0)continue;
1402         if(part->Pt()<fTrackPtCut)continue;
1403         if(kTrackAODMCCharged){
1404           list->Add(part);
1405         }
1406         else {
1407           if(TMath::Abs(part->Eta())>0.9)continue;
1408           list->Add(part);
1409         }
1410         iCount++;
1411       }
1412     }
1413   }// AODMCparticle
1414   list->Sort();
1415   return iCount;
1416 }
1417
1418 /*
1419 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1420   for(int i = 0; i < particles.GetEntries(); i++){
1421     AliVParticle *vp = (AliVParticle*)particles.At(i);
1422     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1423     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1424     jInp.set_user_index(i);
1425     inputParticles.push_back(jInp);
1426   }
1427
1428   return 0;
1429
1430 }
1431 */