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