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 <TRandom.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TRefArray.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetCluster.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliJetFinder.h"
  43 #include "AliJetHeader.h"
  44 #include "AliJetReader.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODTrack.h"
  49 #include "AliAODJet.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliGenPythiaEventHeader.h"
  55 #include "AliJetKineReaderHeader.h"
  56 #include "AliGenCocktailEventHeader.h"
  57 #include "AliInputEventHandler.h"
  58 #include "AliAODJetEventBackground.h"
  59
  60 #include "fastjet/PseudoJet.hh"
  61 #include "fastjet/ClusterSequenceArea.hh"
  62 #include "fastjet/AreaDefinition.hh"
  63 #include "fastjet/JetDefinition.hh"
  64 // get info on how fastjet was configured
  65 #include "fastjet/config.h"
  66
  67
  68 ClassImp(AliAnalysisTaskJetCluster)
  69
  70 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
  71   delete fRef;
  72 }
  73
  74 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
  75   fAOD(0x0),
  76   fAODExtension(0x0),
  77   fRef(new TRefArray),
  78   fUseAODTrackInput(kFALSE),
  79   fUseAODMCInput(kFALSE),
  80   fUseGlobalSelection(kFALSE),
  81   fFilterMask(0),
  82   fTrackTypeRec(kTrackUndef),
  83   fTrackTypeGen(kTrackUndef),
  84   fNSkipLeadingRan(0),
  85   fAvgTrials(1),
  86   fExternalWeight(1),
  87   fRecEtaWindow(0.5),
  88   fTrackPtCut(0.),
  89   fJetOutputMinPt(1),
  90   fNonStdBranch(""),
  91   fNonStdFile(""),
  92   fRparam(1.0),
  93   fAlgorithm(fastjet::kt_algorithm),
  94   fStrategy(fastjet::Best),
  95   fRecombScheme(fastjet::BIpt_scheme),
  96   fAreaType(fastjet::active_area),
  97   fGhostArea(0.01),
  98   fActiveAreaRepeats(1),
  99   fGhostEtamax(1.5),
 100   fh1Xsec(0x0),
 101   fh1Trials(0x0),
 102   fh1PtHard(0x0),
 103   fh1PtHardNoW(0x0),
 104   fh1PtHardTrials(0x0),
 105   fh1NJetsRec(0x0),
 106   fh1NConstRec(0x0),
 107   fh1NConstLeadingRec(0x0),
 108   fh1PtJetsRecIn(0x0),
 109   fh1PtJetsLeadingRecIn(0x0),
 110   fh1PtJetConstRec(0x0),
 111   fh1PtJetConstLeadingRec(0x0),
 112   fh1PtTracksRecIn(0x0),
 113   fh1PtTracksLeadingRecIn(0x0),
 114   fh1NJetsRecRan(0x0),
 115   fh1NConstRecRan(0x0),
 116   fh1PtJetsLeadingRecInRan(0x0),
 117   fh1NConstLeadingRecRan(0x0),
 118   fh1PtJetsRecInRan(0x0),
 119   fh1PtTracksGenIn(0x0),
 120   fh1Nch(0x0),
 121   fh2NRecJetsPt(0x0),
 122   fh2NRecTracksPt(0x0),
 123   fh2NConstPt(0x0),
 124   fh2NConstLeadingPt(0x0),
 125   fh2JetPhiEta(0x0),
 126   fh2LeadingJetPhiEta(0x0),
 127   fh2JetEtaPt(0x0),
 128   fh2LeadingJetEtaPt(0x0),
 129   fh2TrackEtaPt(0x0),
 130   fh2LeadingTrackEtaPt(0x0),
 131   fh2JetsLeadingPhiEta(0x0),
 132   fh2JetsLeadingPhiPt(0x0),
 133   fh2TracksLeadingPhiEta(0x0),
 134   fh2TracksLeadingPhiPt(0x0),
 135   fh2TracksLeadingJetPhiPt(0x0),
 136   fh2JetsLeadingPhiPtW(0x0),
 137   fh2TracksLeadingPhiPtW(0x0),
 138   fh2TracksLeadingJetPhiPtW(0x0),
 139   fh2NRecJetsPtRan(0x0),
 140   fh2NConstPtRan(0x0),
 141   fh2NConstLeadingPtRan(0x0),
 142   fh2PtNch(0x0),
 143   fh2PtNchRan(0x0),
 144   fh2PtNchN(0x0),
 145   fh2PtNchNRan(0x0),
 146   fh2TracksLeadingJetPhiPtRan(0x0),
 147   fh2TracksLeadingJetPhiPtWRan(0x0),
 148   fHistList(0x0)
 149 {
 150
 151 }
 152
 153 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 154   AliAnalysisTaskSE(name),
 155   fAOD(0x0),
 156   fAODExtension(0x0),
 157   fRef(new TRefArray),
 158   fUseAODTrackInput(kFALSE),
 159   fUseAODMCInput(kFALSE),
 160   fUseGlobalSelection(kFALSE),
 161   fFilterMask(0),
 162   fTrackTypeRec(kTrackUndef),
 163   fTrackTypeGen(kTrackUndef),
 164   fNSkipLeadingRan(0),
 165   fAvgTrials(1),
 166   fExternalWeight(1),
 167   fRecEtaWindow(0.5),
 168   fTrackPtCut(0.),
 169   fJetOutputMinPt(1),
 170   fNonStdBranch(""),
 171   fNonStdFile(""),
 172   fRparam(1.0),
 173   fAlgorithm(fastjet::kt_algorithm),
 174   fStrategy(fastjet::Best),
 175   fRecombScheme(fastjet::BIpt_scheme),
 176   fAreaType(fastjet::active_area),
 177   fGhostArea(0.01),
 178   fActiveAreaRepeats(1),
 179   fGhostEtamax(1.5),
 180   fh1Xsec(0x0),
 181   fh1Trials(0x0),
 182   fh1PtHard(0x0),
 183   fh1PtHardNoW(0x0),
 184   fh1PtHardTrials(0x0),
 185   fh1NJetsRec(0x0),
 186   fh1NConstRec(0x0),
 187   fh1NConstLeadingRec(0x0),
 188   fh1PtJetsRecIn(0x0),
 189   fh1PtJetsLeadingRecIn(0x0),
 190   fh1PtJetConstRec(0x0),
 191   fh1PtJetConstLeadingRec(0x0),
 192   fh1PtTracksRecIn(0x0),
 193   fh1PtTracksLeadingRecIn(0x0),
 194   fh1NJetsRecRan(0x0),
 195   fh1NConstRecRan(0x0),
 196   fh1PtJetsLeadingRecInRan(0x0),
 197   fh1NConstLeadingRecRan(0x0),
 198   fh1PtJetsRecInRan(0x0),
 199   fh1PtTracksGenIn(0x0),
 200   fh1Nch(0x0),
 201   fh2NRecJetsPt(0x0),
 202   fh2NRecTracksPt(0x0),
 203   fh2NConstPt(0x0),
 204   fh2NConstLeadingPt(0x0),
 205   fh2JetPhiEta(0x0),
 206   fh2LeadingJetPhiEta(0x0),
 207   fh2JetEtaPt(0x0),
 208   fh2LeadingJetEtaPt(0x0),
 209   fh2TrackEtaPt(0x0),
 210   fh2LeadingTrackEtaPt(0x0),
 211   fh2JetsLeadingPhiEta(0x0),
 212   fh2JetsLeadingPhiPt(0x0),
 213   fh2TracksLeadingPhiEta(0x0),
 214   fh2TracksLeadingPhiPt(0x0),
 215   fh2TracksLeadingJetPhiPt(0x0),
 216   fh2JetsLeadingPhiPtW(0x0),
 217   fh2TracksLeadingPhiPtW(0x0),
 218   fh2TracksLeadingJetPhiPtW(0x0),
 219   fh2NRecJetsPtRan(0x0),
 220   fh2NConstPtRan(0x0),
 221   fh2NConstLeadingPtRan(0x0),
 222   fh2PtNch(0x0),
 223   fh2PtNchRan(0x0),
 224   fh2PtNchN(0x0),
 225   fh2PtNchNRan(0x0),
 226   fh2TracksLeadingJetPhiPtRan(0x0),
 227   fh2TracksLeadingJetPhiPtWRan(0x0),
 228   fHistList(0x0)
 229 {
 230   DefineOutput(1, TList::Class());
 231 }
 232
 233
 234
 235 Bool_t AliAnalysisTaskJetCluster::Notify()
 236 {
 237   //
 238   // Implemented Notify() to read the cross sections
 239   // and number of trials from pyxsec.root
 240   //
 241   return kTRUE;
 242 }
 243
 244 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 245 {
 246
 247   //
 248   // Create the output container
 249   //
 250
 251
 252
 253
 254   // Connect the AOD
 255
 256
 257   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 258
 259
 260
 261   if(fNonStdBranch.Length()!=0)
 262     {
 263       // only create the output branch if we have a name
 264       // Create a new branch for jets...
 265       //  -> cleared in the UserExec....
 266       // here we can also have the case that the brnaches are written to a separate file
 267
 268       TClonesArray *tca = new TClonesArray("AliAODJet", 0);
 269       tca->SetName(fNonStdBranch.Data());
 270       AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
 271
 272
 273       TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
 274       tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 275       AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
 276
 277
 278      if(!AODEvent() || !AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 279         AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
 280         evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 281         AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
 282
 283      }
 284
 285       if(fNonStdFile.Length()!=0){
 286         //
 287         // case that we have an AOD extension we need to fetch the jets from the extended output
 288         // we identifay the extension aod event by looking for the branchname
 289         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 290         TObjArray* extArray = aodH->GetExtensions();
 291         if (extArray) {
 292           TIter next(extArray);
 293           while ((fAODExtension=(AliAODExtension*)next())){
 294             TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
 295             if(fDebug>10){
 296               Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
 297               fAODExtension->GetAOD()->Dump();
 298             }
 299             if(obj){
 300               if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
 301               break;
 302             }
 303             fAODExtension = 0;
 304           }
 305         }
 306       }
 307     }
 308
 309
 310   OpenFile(1);
 311   if(!fHistList)fHistList = new TList();
 312
 313   Bool_t oldStatus = TH1::AddDirectoryStatus();
 314   TH1::AddDirectory(kFALSE);
 315
 316   //
 317   //  Histogram
 318
 319   const Int_t nBinPt = 200;
 320   Double_t binLimitsPt[nBinPt+1];
 321   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 322     if(iPt == 0){
 323       binLimitsPt[iPt] = 0.0;
 324     }
 325     else {// 1.0
 326       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 0.5;
 327     }
 328   }
 329
 330   const Int_t nBinPhi = 90;
 331   Double_t binLimitsPhi[nBinPhi+1];
 332   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 333     if(iPhi==0){
 334       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 335     }
 336     else{
 337       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 338     }
 339   }
 340
 341
 342
 343   const Int_t nBinEta = 40;
 344   Double_t binLimitsEta[nBinEta+1];
 345   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 346     if(iEta==0){
 347       binLimitsEta[iEta] = -2.0;
 348     }
 349     else{
 350       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 351     }
 352   }
 353
 354   const int nChMax = 100;
 355
 356   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 357   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 358
 359   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 360   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 361
 362
 363   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 364   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 365
 366   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 367   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 368   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 369   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 370
 371
 372   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 373   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 374   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 375
 376   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 377   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 378   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 379   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 380   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 381   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 382   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 383   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 384   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 385   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 386
 387   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 388   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 389   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 390   //
 391
 392
 393   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 394   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 395   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 396   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 397
 398   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 399   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);
 400   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);
 401   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);
 402
 403
 404
 405   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 406                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 407   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 408                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 409
 410   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 411                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 412   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 413                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 414
 415   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 416                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 417   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 418                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 419
 420
 421
 422   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 423                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 424   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 425                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 426   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 427                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 428   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 429                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 430   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 431                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 432   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 433                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 434
 435   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 436                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 437   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 438                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 439
 440   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 441                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 442   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 443                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 444
 445
 446
 447   const Int_t saveLevel = 3; // large save level more histos
 448   if(saveLevel>0){
 449     fHistList->Add(fh1Xsec);
 450     fHistList->Add(fh1Trials);
 451
 452     fHistList->Add(fh1NJetsRec);
 453     fHistList->Add(fh1NConstRec);
 454     fHistList->Add(fh1NConstLeadingRec);
 455     fHistList->Add(fh1PtJetsRecIn);
 456     fHistList->Add(fh1PtJetsLeadingRecIn);
 457     fHistList->Add(fh1PtTracksRecIn);
 458     fHistList->Add(fh1PtTracksLeadingRecIn);
 459     fHistList->Add(fh1PtJetConstRec);
 460     fHistList->Add(fh1PtJetConstLeadingRec);
 461     fHistList->Add(fh1NJetsRecRan);
 462     fHistList->Add(fh1NConstRecRan);
 463     fHistList->Add(fh1PtJetsLeadingRecInRan);
 464     fHistList->Add(fh1NConstLeadingRecRan);
 465     fHistList->Add(fh1PtJetsRecInRan);
 466     fHistList->Add(fh1Nch);
 467     fHistList->Add(fh2NRecJetsPt);
 468     fHistList->Add(fh2NRecTracksPt);
 469     fHistList->Add(fh2NConstPt);
 470     fHistList->Add(fh2NConstLeadingPt);
 471     fHistList->Add(fh2PtNch);
 472     fHistList->Add(fh2PtNchRan);
 473     fHistList->Add(fh2PtNchN);
 474     fHistList->Add(fh2PtNchNRan);
 475     fHistList->Add(fh2JetPhiEta);
 476     fHistList->Add(fh2LeadingJetPhiEta);
 477     fHistList->Add(fh2JetEtaPt);
 478     fHistList->Add(fh2LeadingJetEtaPt);
 479     fHistList->Add(fh2TrackEtaPt);
 480     fHistList->Add(fh2LeadingTrackEtaPt);
 481     fHistList->Add(fh2JetsLeadingPhiEta );
 482     fHistList->Add(fh2JetsLeadingPhiPt);
 483     fHistList->Add(fh2TracksLeadingPhiEta);
 484     fHistList->Add(fh2TracksLeadingPhiPt);
 485     fHistList->Add(fh2TracksLeadingJetPhiPt);
 486     fHistList->Add(fh2JetsLeadingPhiPtW);
 487     fHistList->Add(fh2TracksLeadingPhiPtW);
 488     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 489     fHistList->Add(fh2NRecJetsPtRan);
 490     fHistList->Add(fh2NConstPtRan);
 491     fHistList->Add(fh2NConstLeadingPtRan);
 492     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 493     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 494     }
 495
 496   // =========== Switch on Sumw2 for all histos ===========
 497   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 498     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 499     if (h1){
 500       h1->Sumw2();
 501       continue;
 502     }
 503     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 504     if(hn)hn->Sumw2();
 505   }
 506   TH1::AddDirectory(oldStatus);
 507 }
 508
 509 void AliAnalysisTaskJetCluster::Init()
 510 {
 511   //
 512   // Initialization
 513   //
 514
 515   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 516
 517 }
 518
 519 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 520 {
 521
 522   if(fUseGlobalSelection){
 523     // no selection by the service task, we continue
 524     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 525     PostData(1, fHistList);
 526     return;
 527   }
 528
 529   // handle and reset the output jet branch
 530   // only need this once
 531   TClonesArray* jarray = 0;
 532   TClonesArray* jarrayran = 0;
 533   AliAODJetEventBackground*  evBkg = 0;
 534   if(fNonStdBranch.Length()!=0) {
 535     if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
 536     if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
 537     if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
 538     if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
 539     if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
 540     if(jarrayran)jarrayran->Delete();    // this is our responsibility, clear before filling again
 541
 542     if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
 543     if(!evBkg)  evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
 544     if(evBkg)evBkg->Reset();
 545
 546   }
 547
 548
 549
 550   //
 551   // Execute analysis for current event
 552   //
 553   AliESDEvent *fESD = 0;
 554   if(fUseAODTrackInput){
 555     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 556     if(!fAOD){
 557       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 558       return;
 559     }
 560     // fethc the header
 561   }
 562   else{
 563     //  assume that the AOD is in the general output...
 564     fAOD  = AODEvent();
 565     if(!fAOD){
 566       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 567       return;
 568     }
 569     if(fDebug>0){
 570       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 571     }
 572   }
 573
 574   Bool_t selectEvent =  false;
 575   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 576
 577   if(fAOD){
 578     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 579     TString vtxTitle(vtxAOD->GetTitle());
 580
 581     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 582         Float_t zvtx = vtxAOD->GetZ();
 583         Float_t yvtx = vtxAOD->GetY();
 584         Float_t xvtx = vtxAOD->GetX();
 585         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 586         if(TMath::Abs(zvtx)<8.&&r2<1.){
 587           if(physicsSelection)selectEvent = true;
 588         }
 589     }
 590   }
 591   if(!selectEvent){
 592     PostData(1, fHistList);
 593     return;
 594   }
 595
 596   fh1Trials->Fill("#sum{ntrials}",1);
 597
 598
 599   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 600
 601   // ==== General variables needed
 602
 603
 604
 605   // we simply fetch the tracks/mc particles as a list of AliVParticles
 606
 607   TList recParticles;
 608   TList genParticles;
 609
 610   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 611   Float_t nCh = recParticles.GetEntries();
 612   fh1Nch->Fill(nCh);
 613   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 614   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 615   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 616
 617   // find the jets....
 618
 619   vector<fastjet::PseudoJet> inputParticlesRec;
 620   vector<fastjet::PseudoJet> inputParticlesRecRan;
 621   for(int i = 0; i < recParticles.GetEntries(); i++){
 622     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 623     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 624     // we talk total momentum here
 625     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 626     jInp.set_user_index(i);
 627     inputParticlesRec.push_back(jInp);
 628
 629     // the randomized input changes eta and phi, but keeps the p_T
 630     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 631       Double_t pT = vp->Pt();
 632       Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
 633       Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
 634
 635       Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
 636       Double_t pZ = pT/TMath::Tan(theta);
 637
 638       Double_t pX = pT * TMath::Cos(phi);
 639       Double_t pY = pT * TMath::Sin(phi);
 640       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 641       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 642       jInpRan.set_user_index(i);
 643       inputParticlesRecRan.push_back(jInpRan);
 644     }
 645
 646     // fill the tref array, only needed when we write out jets
 647     if(jarray){
 648       if(i == 0){
 649         fRef->Delete(); // make sure to delete before placement new...
 650         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
 651       }
 652       fRef->Add(vp);
 653     }
 654   }
 655
 656   if(inputParticlesRec.size()==0){
 657     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 658     PostData(1, fHistList);
 659     return;
 660   }
 661
 662   // run fast jet
 663   // employ setters for these...
 664
 665
 666   // now create the object that holds info about ghosts
 667   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 668   fastjet::AreaType areaType =   fastjet::active_area;
 669   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 670     fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 671   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 672   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 673
 674   //range where to compute background
 675   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 676   phiMin = 0;
 677   phiMax = 2*TMath::Pi();
 678   rapMax = fGhostEtamax - fRparam;
 679   rapMin = - fGhostEtamax + fRparam;
 680   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 681
 682
 683
 684   inclusiveJets = clustSeq.inclusive_jets();
 685   sortedJets    = sorted_by_pt(inclusiveJets);
 686
 687   fh1NJetsRec->Fill(sortedJets.size());
 688
 689  // loop over all jets an fill information, first on is the leading jet
 690
 691   Int_t nRecOver = inclusiveJets.size();
 692   Int_t nRec     = inclusiveJets.size();
 693   if(inclusiveJets.size()>0){
 694     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 695     Float_t pt = leadingJet.Pt();
 696     Int_t nAodOutJets = 0;
 697     Int_t nAodOutTracks = 0;
 698     AliAODJet *aodOutJet = 0;
 699
 700     Int_t iCount = 0;
 701     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 702       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 703       while(pt<ptCut&&iCount<nRec){
 704         nRecOver--;
 705         iCount++;
 706         if(iCount<nRec){
 707           pt = sortedJets[iCount].perp();
 708         }
 709       }
 710       if(nRecOver<=0)break;
 711       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 712     }
 713     Float_t phi = leadingJet.Phi();
 714     if(phi<0)phi+=TMath::Pi()*2.;
 715     Float_t eta = leadingJet.Eta();
 716     pt = leadingJet.Pt();
 717
 718     // correlation of leading jet with tracks
 719     TIterator *recIter = recParticles.MakeIterator();
 720     recIter->Reset();
 721     AliVParticle *tmpRecTrack = 0;
 722     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 723       Float_t tmpPt = tmpRecTrack->Pt();
 724       // correlation
 725       Float_t tmpPhi =  tmpRecTrack->Phi();
 726       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 727       Float_t dPhi = phi - tmpPhi;
 728       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 729       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 730       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 731       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 732     }
 733
 734
 735
 736    for(int j = 0; j < nRec;j++){
 737      AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 738      aodOutJet = 0;
 739      nAodOutTracks = 0;
 740      Float_t tmpPt = tmpRec.Pt();
 741      fh1PtJetsRecIn->Fill(tmpPt);
 742      // Fill Spectra with constituents
 743      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 744      fh1NConstRec->Fill(constituents.size());
 745      fh2PtNch->Fill(nCh,tmpPt);
 746      fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 747      fh2NConstPt->Fill(tmpPt,constituents.size());
 748      // loop over constiutents and fill spectrum
 749
 750      // Add the jet information and the track references to the output AOD
 751
 752      if(tmpPt>fJetOutputMinPt&&jarray){
 753        aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
 754        Double_t area=clustSeq.area(sortedJets[j]);
 755
 756        aodOutJet->SetEffArea(area,0);
 757      }
 758
 759
 760      for(unsigned int ic = 0; ic < constituents.size();ic++){
 761        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 762        fh1PtJetConstRec->Fill(part->Pt());
 763        if(aodOutJet){
 764          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
 765        }
 766        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
 767      }
 768
 769
 770
 771
 772
 773
 774
 775
 776      // correlation
 777      Float_t tmpPhi =  tmpRec.Phi();
 778      Float_t tmpEta =  tmpRec.Eta();
 779      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 780
 781      if(j==0){
 782        fh1PtJetsLeadingRecIn->Fill(tmpPt);
 783        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 784        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 785        fh1NConstLeadingRec->Fill(constituents.size());
 786        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 787        continue;
 788      }
 789      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
 790      fh2JetEtaPt->Fill(tmpEta,tmpPt);
 791      Float_t dPhi = phi - tmpPhi;
 792      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 793      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 794      Float_t dEta = eta - tmpRec.Eta();
 795      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 796      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 797      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 798    }
 799    delete recIter;
 800
 801       //background estimates:all bckg jets(0) & wo the 2 hardest(1)
 802
 803    if(evBkg){
 804      vector<fastjet::PseudoJet> jets2=sortedJets;
 805      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
 806      Double_t bkg1=0;
 807      Double_t sigma1=0.;
 808      Double_t meanarea1=0.;
 809      Double_t bkg2=0;
 810      Double_t sigma2=0.;
 811      Double_t meanarea2=0.;
 812
 813      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
 814      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
 815      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
 816      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
 817    }
 818
 819
 820
 821   }
 822
 823  // fill track information
 824  Int_t nTrackOver = recParticles.GetSize();
 825   // do the same for tracks and jets
 826  if(nTrackOver>0){
 827    TIterator *recIter = recParticles.MakeIterator();
 828    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 829    Float_t pt = tmpRec->Pt();
 830    //    Printf("Leading track p_t %3.3E",pt);
 831    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 832      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 833      while(pt<ptCut&&tmpRec){
 834        nTrackOver--;
 835        tmpRec = (AliVParticle*)(recIter->Next());
 836        if(tmpRec){
 837          pt = tmpRec->Pt();
 838        }
 839      }
 840      if(nTrackOver<=0)break;
 841      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 842    }
 843
 844    recIter->Reset();
 845    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 846    Float_t phi = leading->Phi();
 847    if(phi<0)phi+=TMath::Pi()*2.;
 848    Float_t eta = leading->Eta();
 849    pt  = leading->Pt();
 850    while((tmpRec = (AliVParticle*)(recIter->Next()))){
 851      Float_t tmpPt = tmpRec->Pt();
 852      Float_t tmpEta = tmpRec->Eta();
 853      fh1PtTracksRecIn->Fill(tmpPt);
 854      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
 855      if(tmpRec==leading){
 856        fh1PtTracksLeadingRecIn->Fill(tmpPt);
 857        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
 858        continue;
 859      }
 860       // correlation
 861      Float_t tmpPhi =  tmpRec->Phi();
 862
 863      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 864      Float_t dPhi = phi - tmpPhi;
 865      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 866      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 867      Float_t dEta = eta - tmpRec->Eta();
 868      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 869      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 870      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 871    }
 872    delete recIter;
 873  }
 874
 875  // find the random jets
 876  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
 877  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
 878
 879  inclusiveJetsRan = clustSeqRan.inclusive_jets();
 880  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
 881
 882  // fill the jet information from random track
 883
 884   fh1NJetsRecRan->Fill(sortedJetsRan.size());
 885
 886  // loop over all jets an fill information, first on is the leading jet
 887
 888  Int_t nRecOverRan = inclusiveJetsRan.size();
 889  Int_t nRecRan     = inclusiveJetsRan.size();
 890  if(inclusiveJetsRan.size()>0){
 891    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
 892    Float_t pt = leadingJet.Pt();
 893
 894    Int_t iCount = 0;
 895    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
 896      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
 897       while(pt<ptCut&&iCount<nRecRan){
 898         nRecOverRan--;
 899         iCount++;
 900         if(iCount<nRecRan){
 901           pt = sortedJetsRan[iCount].perp();
 902         }
 903       }
 904       if(nRecOverRan<=0)break;
 905       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
 906     }
 907     Float_t phi = leadingJet.Phi();
 908     if(phi<0)phi+=TMath::Pi()*2.;
 909     pt = leadingJet.Pt();
 910
 911     // correlation of leading jet with random tracks
 912
 913     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
 914       {
 915         Float_t tmpPt = inputParticlesRecRan[ip].perp();
 916         // correlation
 917         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
 918         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 919         Float_t dPhi = phi - tmpPhi;
 920         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 921         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 922         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
 923         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
 924       }
 925
 926     Int_t nAodOutJetsRan = 0;
 927      AliAODJet *aodOutJetRan = 0;
 928     for(int j = 0; j < nRecRan;j++){
 929       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
 930       Float_t tmpPt = tmpRec.Pt();
 931       fh1PtJetsRecInRan->Fill(tmpPt);
 932       // Fill Spectra with constituents
 933       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
 934       fh1NConstRecRan->Fill(constituents.size());
 935       fh2NConstPtRan->Fill(tmpPt,constituents.size());
 936       fh2PtNchRan->Fill(nCh,tmpPt);
 937       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
 938
 939
 940      if(tmpPt>fJetOutputMinPt&&jarrayran){
 941        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
 942        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
 943
 944        aodOutJetRan->SetEffArea(arearan,0);    }
 945
 946
 947
 948
 949      for(unsigned int ic = 0; ic < constituents.size();ic++){
 950        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 951        // fh1PtJetConstRec->Fill(part->Pt());
 952        if(aodOutJetRan){
 953          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
 954        }}
 955
 956
 957
 958
 959       // correlation
 960       Float_t tmpPhi =  tmpRec.Phi();
 961       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 962
 963       if(j==0){
 964         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
 965         fh1NConstLeadingRecRan->Fill(constituents.size());
 966         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
 967         continue;
 968       }
 969     }
 970
 971
 972     if(evBkg){
 973      Double_t bkg3=0.;
 974      Double_t sigma3=0.;
 975      Double_t meanarea3=0.;
 976      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
 977      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
 978     }
 979
 980
 981
 982  }
 983
 984
 985  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 986  PostData(1, fHistList);
 987 }
 988
 989 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
 990 {
 991 // Terminate analysis
 992 //
 993     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 994 }
 995
 996
 997 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 998
 999   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1000
1001   Int_t iCount = 0;
1002   if(type==kTrackAOD){
1003     AliAODEvent *aod = 0;
1004     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1005     else aod = AODEvent();
1006     if(!aod){
1007       return iCount;
1008     }
1009     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1010       AliAODTrack *tr = aod->GetTrack(it);
1011       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1012       if(TMath::Abs(tr->Eta())>0.9)continue;
1013       //      if(tr->Pt()<0.3)continue;
1014       list->Add(tr);
1015       iCount++;
1016     }
1017   }
1018   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1019     AliMCEvent* mcEvent = MCEvent();
1020     if(!mcEvent)return iCount;
1021     // we want to have alivpartilces so use get track
1022     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1023       if(!mcEvent->IsPhysicalPrimary(it))continue;
1024       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1025       if(type == kTrackKineAll){
1026         list->Add(part);
1027         iCount++;
1028       }
1029       else if(type == kTrackKineCharged){
1030         if(part->Particle()->GetPDG()->Charge()==0)continue;
1031         list->Add(part);
1032         iCount++;
1033       }
1034     }
1035   }
1036   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1037     AliAODEvent *aod = 0;
1038     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1039     else aod = AODEvent();
1040     if(!aod)return iCount;
1041     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1042     if(!tca)return iCount;
1043     for(int it = 0;it < tca->GetEntriesFast();++it){
1044       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1045       if(!part->IsPhysicalPrimary())continue;
1046       if(type == kTrackAODMCAll){
1047         list->Add(part);
1048         iCount++;
1049       }
1050       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1051         if(part->Charge()==0)continue;
1052         if(kTrackAODMCCharged){
1053           list->Add(part);
1054         }
1055         else {
1056           if(TMath::Abs(part->Eta())>0.9)continue;
1057           list->Add(part);
1058         }
1059         iCount++;
1060       }
1061     }
1062   }// AODMCparticle
1063   list->Sort();
1064   return iCount;
1065
1066 }
1067
1068 /*
1069 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1070   for(int i = 0; i < particles.GetEntries(); i++){
1071     AliVParticle *vp = (AliVParticle*)particles.At(i);
1072     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1073     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1074     jInp.set_user_index(i);
1075     inputParticles.push_back(jInp);
1076   }
1077
1078   return 0;
1079
1080 }
1081 */