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      Float_t tmpPtBack = 0;
 969      if(externalBackground){
 970        // carefull has to be filled in a task before
 971        // todo, ReArrange to the botom
 972        tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
 973      }
 974      tmpPt = tmpPt - tmpPtBack;
 975      if(tmpPt<0)tmpPt = 0; // avoid negative weights...
 976
 977      fh1PtJetsRecIn->Fill(tmpPt);
 978      // Fill Spectra with constituents
 979      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 980      fh1NConstRec->Fill(constituents.size());
 981      fh2PtNch->Fill(nCh,tmpPt);
 982      fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 983      fh2NConstPt->Fill(tmpPt,constituents.size());
 984      // loop over constiutents and fill spectrum
 985
 986      // Add the jet information and the track references to the output AOD
 987
 988      if(tmpPt>fJetOutputMinPt&&jarray){
 989        aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
 990        Double_t area1 = clustSeq.area(sortedJets[j]);
 991        aodOutJet->SetEffArea(area1,0);
 992      }
 993
 994
 995      for(unsigned int ic = 0; ic < constituents.size();ic++){
 996        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 997        fh1PtJetConstRec->Fill(part->Pt());
 998        if(aodOutJet){
 999          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1000        }
1001        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1002      }
1003
1004      // correlation
1005      Float_t tmpPhi =  tmpRec.Phi();
1006      Float_t tmpEta =  tmpRec.Eta();
1007      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1008
1009
1010
1011      if(j==0){
1012        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1013        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1014        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1015        fh1NConstLeadingRec->Fill(constituents.size());
1016        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1017        continue;
1018      }
1019      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1020      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1021      Float_t dPhi = phi - tmpPhi;
1022      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1023      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1024      Float_t dEta = eta - tmpRec.Eta();
1025      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1026      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1027      if(cenClass>=0){
1028        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1029        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1030      }
1031      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1032    }
1033    delete recIter;
1034
1035    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1036
1037    if(evBkg){
1038      vector<fastjet::PseudoJet> jets2=sortedJets;
1039      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1040      Double_t bkg1=0;
1041      Double_t sigma1=0.;
1042      Double_t meanarea1=0.;
1043      Double_t bkg2=0;
1044      Double_t sigma2=0.;
1045      Double_t meanarea2=0.;
1046
1047      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
1048      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1049
1050      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1051      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1052
1053      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1054      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1055      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1056      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1057
1058    }
1059   }
1060
1061
1062
1063
1064
1065   // fill track information
1066   Int_t nTrackOver = recParticles.GetSize();
1067   // do the same for tracks and jets
1068
1069   if(nTrackOver>0){
1070    TIterator *recIter = recParticles.MakeIterator();
1071    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1072    Float_t pt = tmpRec->Pt();
1073
1074    //    Printf("Leading track p_t %3.3E",pt);
1075    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1076      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1077      while(pt<ptCut&&tmpRec){
1078        nTrackOver--;
1079        tmpRec = (AliVParticle*)(recIter->Next());
1080        if(tmpRec){
1081          pt = tmpRec->Pt();
1082        }
1083      }
1084      if(nTrackOver<=0)break;
1085      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1086    }
1087
1088    recIter->Reset();
1089    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1090    Float_t phi = leading->Phi();
1091    if(phi<0)phi+=TMath::Pi()*2.;
1092    Float_t eta = leading->Eta();
1093    pt  = leading->Pt();
1094    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1095      Float_t tmpPt = tmpRec->Pt();
1096      Float_t tmpEta = tmpRec->Eta();
1097      fh1PtTracksRecIn->Fill(tmpPt);
1098      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1099      if(tmpRec==leading){
1100        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1101        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1102        continue;
1103      }
1104       // correlation
1105      Float_t tmpPhi =  tmpRec->Phi();
1106
1107      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1108      Float_t dPhi = phi - tmpPhi;
1109      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1110      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1111      Float_t dEta = eta - tmpRec->Eta();
1112      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1113      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1114      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1115    }
1116    delete recIter;
1117  }
1118
1119  // find the random jets
1120  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1121  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1122
1123  inclusiveJetsRan = clustSeqRan.inclusive_jets();
1124  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1125
1126  // fill the jet information from random track
1127
1128   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1129
1130  // loop over all jets an fill information, first on is the leading jet
1131
1132  Int_t nRecOverRan = inclusiveJetsRan.size();
1133  Int_t nRecRan     = inclusiveJetsRan.size();
1134  if(inclusiveJetsRan.size()>0){
1135    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1136    Float_t pt = leadingJet.Pt();
1137
1138    Int_t iCount = 0;
1139    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1140      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1141       while(pt<ptCut&&iCount<nRecRan){
1142         nRecOverRan--;
1143         iCount++;
1144         if(iCount<nRecRan){
1145           pt = sortedJetsRan[iCount].perp();
1146         }
1147       }
1148       if(nRecOverRan<=0)break;
1149       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1150     }
1151     Float_t phi = leadingJet.Phi();
1152     if(phi<0)phi+=TMath::Pi()*2.;
1153     pt = leadingJet.Pt();
1154
1155     // correlation of leading jet with random tracks
1156
1157     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1158       {
1159         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1160         // correlation
1161         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1162         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1163         Float_t dPhi = phi - tmpPhi;
1164         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1165         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1166         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1167         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1168       }
1169
1170     Int_t nAodOutJetsRan = 0;
1171      AliAODJet *aodOutJetRan = 0;
1172     for(int j = 0; j < nRecRan;j++){
1173       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1174       Float_t tmpPt = tmpRec.Pt();
1175       fh1PtJetsRecInRan->Fill(tmpPt);
1176       // Fill Spectra with constituents
1177       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1178       fh1NConstRecRan->Fill(constituents.size());
1179       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1180       fh2PtNchRan->Fill(nCh,tmpPt);
1181       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1182
1183
1184      if(tmpPt>fJetOutputMinPt&&jarrayran){
1185        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1186        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1187
1188        aodOutJetRan->SetEffArea(arearan,0);    }
1189
1190
1191
1192
1193      for(unsigned int ic = 0; ic < constituents.size();ic++){
1194        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1195        // fh1PtJetConstRec->Fill(part->Pt());
1196        if(aodOutJetRan){
1197          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1198        }
1199      }
1200
1201
1202
1203
1204       // correlation
1205       Float_t tmpPhi =  tmpRec.Phi();
1206       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1207
1208       if(j==0){
1209         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1210         fh1NConstLeadingRecRan->Fill(constituents.size());
1211         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1212         continue;
1213       }
1214     }
1215
1216
1217     if(evBkg){
1218      Double_t bkg3=0.;
1219      Double_t sigma3=0.;
1220      Double_t meanarea3=0.;
1221      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1222      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1223      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1224      /*
1225      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1226      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1227      */
1228     }
1229
1230
1231
1232  }
1233
1234
1235  // do the event selection if activated
1236  if(fJetTriggerPtCut>0){
1237    bool select = false;
1238    Float_t minPt = fJetTriggerPtCut;
1239    /*
1240    // hard coded for now ...
1241    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1242    if(cent<10)minPt = 50;
1243    else if(cent<30)minPt = 42;
1244    else if(cent<50)minPt = 28;
1245    else if(cent<80)minPt = 18;
1246    */
1247    float rho = 0;
1248    if(externalBackground)rho = externalBackground->GetBackground(2);
1249    if(jarray){
1250      for(int i = 0;i < jarray->GetEntriesFast();i++){
1251        AliAODJet *jet = (AliAODJet*)jarray->At(i);
1252        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1253        if(ptSub>=minPt){
1254          select = true;
1255          break;
1256        }
1257      }
1258    }
1259
1260    if(select){
1261      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1262      fh1CentralitySelect->Fill(cent);
1263      aodH->SetFillAOD(kTRUE);
1264    }
1265  }
1266
1267  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1268  PostData(1, fHistList);
1269 }
1270
1271 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1272 {
1273 // Terminate analysis
1274 //
1275     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1276 }
1277
1278
1279 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1280
1281   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1282
1283   Int_t iCount = 0;
1284   if(type==kTrackAOD){
1285     AliAODEvent *aod = 0;
1286     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1287     else aod = AODEvent();
1288     if(!aod){
1289       return iCount;
1290     }
1291     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1292       AliAODTrack *tr = aod->GetTrack(it);
1293       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1294       if(TMath::Abs(tr->Eta())>0.9)continue;
1295       if(tr->Pt()<fTrackPtCut)continue;
1296       list->Add(tr);
1297       iCount++;
1298     }
1299   }
1300   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1301     AliMCEvent* mcEvent = MCEvent();
1302     if(!mcEvent)return iCount;
1303     // we want to have alivpartilces so use get track
1304     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1305       if(!mcEvent->IsPhysicalPrimary(it))continue;
1306       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1307       if(type == kTrackKineAll){
1308         if(part->Pt()<fTrackPtCut)continue;
1309         list->Add(part);
1310         iCount++;
1311       }
1312       else if(type == kTrackKineCharged){
1313         if(part->Particle()->GetPDG()->Charge()==0)continue;
1314         if(part->Pt()<fTrackPtCut)continue;
1315         list->Add(part);
1316         iCount++;
1317       }
1318     }
1319   }
1320   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1321     AliAODEvent *aod = 0;
1322     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1323     else aod = AODEvent();
1324     if(!aod)return iCount;
1325     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1326     if(!tca)return iCount;
1327     for(int it = 0;it < tca->GetEntriesFast();++it){
1328       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1329       if(!part->IsPhysicalPrimary())continue;
1330       if(type == kTrackAODMCAll){
1331         if(part->Pt()<fTrackPtCut)continue;
1332         list->Add(part);
1333         iCount++;
1334       }
1335       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1336         if(part->Charge()==0)continue;
1337         if(part->Pt()<fTrackPtCut)continue;
1338         if(kTrackAODMCCharged){
1339           list->Add(part);
1340         }
1341         else {
1342           if(TMath::Abs(part->Eta())>0.9)continue;
1343           list->Add(part);
1344         }
1345         iCount++;
1346       }
1347     }
1348   }// AODMCparticle
1349   list->Sort();
1350   return iCount;
1351 }
1352
1353 /*
1354 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1355   for(int i = 0; i < particles.GetEntries(); i++){
1356     AliVParticle *vp = (AliVParticle*)particles.At(i);
1357     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1358     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1359     jInp.set_user_index(i);
1360     inputParticles.push_back(jInp);
1361   }
1362
1363   return 0;
1364
1365 }
1366 */