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