JETAN/AliAnalysisTaskJetCluster.cxx

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