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