JETAN/AliAnalysisTaskJetCluster.cxx

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