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