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