JETAN/AliAnalysisTaskJetCluster.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TRandom.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 "AliJetReaderHeader.h"
  46 #include "AliUA1JetHeaderV1.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAODEvent.h"
  49 #include "AliAODHandler.h"
  50 #include "AliAODTrack.h"
  51 #include "AliAODJet.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliJetKineReaderHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliInputEventHandler.h"
  60
  61
  62 #include "fastjet/PseudoJet.hh"
  63 #include "fastjet/ClusterSequenceArea.hh"
  64 #include "fastjet/AreaDefinition.hh"
  65 #include "fastjet/JetDefinition.hh"
  66 // get info on how fastjet was configured
  67 #include "fastjet/config.h"
  68
  69
  70 ClassImp(AliAnalysisTaskJetCluster)
  71
  72 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
  73   fAOD(0x0),
  74   fUseAODTrackInput(kFALSE),
  75   fUseAODMCInput(kFALSE),
  76   fUseGlobalSelection(kFALSE),
  77   fFilterMask(0),
  78   fTrackTypeRec(kTrackUndef),
  79   fTrackTypeGen(kTrackUndef),
  80   fAvgTrials(1),
  81   fExternalWeight(1),
  82   fRecEtaWindow(0.5),
  83   fRparam(1.0),
  84   fAlgorithm(fastjet::kt_algorithm),
  85   fStrategy(fastjet::Best),
  86   fRecombScheme(fastjet::BIpt_scheme),
  87   fAreaType(fastjet::active_area),
  88   fh1Xsec(0x0),
  89   fh1Trials(0x0),
  90   fh1PtHard(0x0),
  91   fh1PtHardNoW(0x0),
  92   fh1PtHardTrials(0x0),
  93   fh1NJetsRec(0x0),
  94   fh1NConstRec(0x0),
  95   fh1NConstLeadingRec(0x0),
  96   fh1PtJetsRecIn(0x0),
  97   fh1PtJetsLeadingRecIn(0x0),
  98   fh1PtJetConstRec(0x0),
  99   fh1PtJetConstLeadingRec(0x0),
 100   fh1PtTracksRecIn(0x0),
 101   fh1PtTracksLeadingRecIn(0x0),
 102   fh1NJetsRecRan(0x0),
 103   fh1NConstRecRan(0x0),
 104   fh1PtJetsLeadingRecInRan(0x0),
 105   fh1NConstLeadingRecRan(0x0),
 106   fh1PtJetsRecInRan(0x0),
 107   fh1PtTracksGenIn(0x0),
 108   fh1Nch(0x0),
 109   fh2NRecJetsPt(0x0),
 110   fh2NRecTracksPt(0x0),
 111   fh2NConstPt(0x0),
 112   fh2NConstLeadingPt(0x0),
 113   fh2JetPhiEta(0x0),
 114   fh2LeadingJetPhiEta(0x0),
 115   fh2JetEtaPt(0x0),
 116   fh2LeadingJetEtaPt(0x0),
 117   fh2TrackEtaPt(0x0),
 118   fh2LeadingTrackEtaPt(0x0),
 119   fh2JetsLeadingPhiEta(0x0),
 120   fh2JetsLeadingPhiPt(0x0),
 121   fh2TracksLeadingPhiEta(0x0),
 122   fh2TracksLeadingPhiPt(0x0),
 123   fh2TracksLeadingJetPhiPt(0x0),
 124   fh2JetsLeadingPhiPtW(0x0),
 125   fh2TracksLeadingPhiPtW(0x0),
 126   fh2TracksLeadingJetPhiPtW(0x0),
 127   fh2NRecJetsPtRan(0x0),
 128   fh2NConstPtRan(0x0),
 129   fh2NConstLeadingPtRan(0x0),
 130   fh2PtNch(0x0),
 131   fh2PtNchRan(0x0),
 132   fh2PtNchN(0x0),
 133   fh2PtNchNRan(0x0),
 134   fh2TracksLeadingJetPhiPtRan(0x0),
 135   fh2TracksLeadingJetPhiPtWRan(0x0),
 136   fHistList(0x0)
 137 {
 138
 139 }
 140
 141 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 142   AliAnalysisTaskSE(name),
 143   fAOD(0x0),
 144   fUseAODTrackInput(kFALSE),
 145   fUseAODMCInput(kFALSE),
 146   fUseGlobalSelection(kFALSE),
 147   fFilterMask(0),
 148   fTrackTypeRec(kTrackUndef),
 149   fTrackTypeGen(kTrackUndef),
 150   fAvgTrials(1),
 151   fExternalWeight(1),
 152   fRecEtaWindow(0.5),
 153   fRparam(1.0),
 154   fAlgorithm(fastjet::kt_algorithm),
 155   fStrategy(fastjet::Best),
 156   fRecombScheme(fastjet::BIpt_scheme),
 157   fAreaType(fastjet::active_area),
 158   fh1Xsec(0x0),
 159   fh1Trials(0x0),
 160   fh1PtHard(0x0),
 161   fh1PtHardNoW(0x0),
 162   fh1PtHardTrials(0x0),
 163   fh1NJetsRec(0x0),
 164   fh1NConstRec(0x0),
 165   fh1NConstLeadingRec(0x0),
 166   fh1PtJetsRecIn(0x0),
 167   fh1PtJetsLeadingRecIn(0x0),
 168   fh1PtJetConstRec(0x0),
 169   fh1PtJetConstLeadingRec(0x0),
 170   fh1PtTracksRecIn(0x0),
 171   fh1PtTracksLeadingRecIn(0x0),
 172   fh1NJetsRecRan(0x0),
 173   fh1NConstRecRan(0x0),
 174   fh1PtJetsLeadingRecInRan(0x0),
 175   fh1NConstLeadingRecRan(0x0),
 176   fh1PtJetsRecInRan(0x0),
 177   fh1PtTracksGenIn(0x0),
 178   fh1Nch(0x0),
 179   fh2NRecJetsPt(0x0),
 180   fh2NRecTracksPt(0x0),
 181   fh2NConstPt(0x0),
 182   fh2NConstLeadingPt(0x0),
 183   fh2JetPhiEta(0x0),
 184   fh2LeadingJetPhiEta(0x0),
 185   fh2JetEtaPt(0x0),
 186   fh2LeadingJetEtaPt(0x0),
 187   fh2TrackEtaPt(0x0),
 188   fh2LeadingTrackEtaPt(0x0),
 189   fh2JetsLeadingPhiEta(0x0),
 190   fh2JetsLeadingPhiPt(0x0),
 191   fh2TracksLeadingPhiEta(0x0),
 192   fh2TracksLeadingPhiPt(0x0),
 193   fh2TracksLeadingJetPhiPt(0x0),
 194   fh2JetsLeadingPhiPtW(0x0),
 195   fh2TracksLeadingPhiPtW(0x0),
 196   fh2TracksLeadingJetPhiPtW(0x0),
 197   fh2NRecJetsPtRan(0x0),
 198   fh2NConstPtRan(0x0),
 199   fh2NConstLeadingPtRan(0x0),
 200   fh2PtNch(0x0),
 201   fh2PtNchRan(0x0),
 202   fh2PtNchN(0x0),
 203   fh2PtNchNRan(0x0),
 204   fh2TracksLeadingJetPhiPtRan(0x0),
 205   fh2TracksLeadingJetPhiPtWRan(0x0),
 206   fHistList(0x0)
 207 {
 208   DefineOutput(1, TList::Class());
 209 }
 210
 211
 212
 213 Bool_t AliAnalysisTaskJetCluster::Notify()
 214 {
 215   //
 216   // Implemented Notify() to read the cross sections
 217   // and number of trials from pyxsec.root
 218   //
 219   return kTRUE;
 220 }
 221
 222 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 223 {
 224
 225   //
 226   // Create the output container
 227   //
 228
 229
 230   // Connect the AOD
 231
 232
 233   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 234
 235   OpenFile(1);
 236   if(!fHistList)fHistList = new TList();
 237
 238   Bool_t oldStatus = TH1::AddDirectoryStatus();
 239   TH1::AddDirectory(kFALSE);
 240
 241   //
 242   //  Histogram
 243
 244   const Int_t nBinPt = 200;
 245   Double_t binLimitsPt[nBinPt+1];
 246   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 247     if(iPt == 0){
 248       binLimitsPt[iPt] = 0.0;
 249     }
 250     else {// 1.0
 251       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 0.5;
 252     }
 253   }
 254
 255   const Int_t nBinPhi = 90;
 256   Double_t binLimitsPhi[nBinPhi+1];
 257   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 258     if(iPhi==0){
 259       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 260     }
 261     else{
 262       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 263     }
 264   }
 265
 266
 267
 268   const Int_t nBinEta = 40;
 269   Double_t binLimitsEta[nBinEta+1];
 270   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 271     if(iEta==0){
 272       binLimitsEta[iEta] = -2.0;
 273     }
 274     else{
 275       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 276     }
 277   }
 278
 279   const int nChMax = 100;
 280
 281   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 282   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 283
 284   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 285   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 286
 287
 288   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 289   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 290
 291   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 292   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 293   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 294   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 295
 296
 297   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 298   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 299   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 300
 301   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 302   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 303   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 304   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 305   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 306   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 307   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 308   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 309   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 310   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 311
 312   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 313   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 314   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 315   //
 316
 317
 318   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 319   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 320   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 321   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 322
 323   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 324   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);
 325   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);
 326   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);
 327
 328
 329
 330   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 331                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 332   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 333                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 334
 335   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 336                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 337   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 338                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 339
 340   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 341                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 342   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 343                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 344
 345
 346
 347   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 348                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 349   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 350                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 351   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 352                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 353   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 354                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 355   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 356                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 357   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 358                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 359
 360   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 361                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 362   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 363                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 364
 365   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 366                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 367   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 368                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 369
 370
 371
 372   const Int_t saveLevel = 3; // large save level more histos
 373   if(saveLevel>0){
 374     fHistList->Add(fh1Xsec);
 375     fHistList->Add(fh1Trials);
 376
 377     fHistList->Add(fh1NJetsRec);
 378     fHistList->Add(fh1NConstRec);
 379     fHistList->Add(fh1NConstLeadingRec);
 380     fHistList->Add(fh1PtJetsRecIn);
 381     fHistList->Add(fh1PtJetsLeadingRecIn);
 382     fHistList->Add(fh1PtTracksRecIn);
 383     fHistList->Add(fh1PtTracksLeadingRecIn);
 384     fHistList->Add(fh1PtJetConstRec);
 385     fHistList->Add(fh1PtJetConstLeadingRec);
 386     fHistList->Add(fh1NJetsRecRan);
 387     fHistList->Add(fh1NConstRecRan);
 388     fHistList->Add(fh1PtJetsLeadingRecInRan);
 389     fHistList->Add(fh1NConstLeadingRecRan);
 390     fHistList->Add(fh1PtJetsRecInRan);
 391     fHistList->Add(fh1Nch);
 392     fHistList->Add(fh2NRecJetsPt);
 393     fHistList->Add(fh2NRecTracksPt);
 394     fHistList->Add(fh2NConstPt);
 395     fHistList->Add(fh2NConstLeadingPt);
 396     fHistList->Add(fh2PtNch);
 397     fHistList->Add(fh2PtNchRan);
 398     fHistList->Add(fh2PtNchN);
 399     fHistList->Add(fh2PtNchNRan);
 400     fHistList->Add(fh2JetPhiEta);
 401     fHistList->Add(fh2LeadingJetPhiEta);
 402     fHistList->Add(fh2JetEtaPt);
 403     fHistList->Add(fh2LeadingJetEtaPt);
 404     fHistList->Add(fh2TrackEtaPt);
 405     fHistList->Add(fh2LeadingTrackEtaPt);
 406     fHistList->Add(fh2JetsLeadingPhiEta );
 407     fHistList->Add(fh2JetsLeadingPhiPt);
 408     fHistList->Add(fh2TracksLeadingPhiEta);
 409     fHistList->Add(fh2TracksLeadingPhiPt);
 410     fHistList->Add(fh2TracksLeadingJetPhiPt);
 411     fHistList->Add(fh2JetsLeadingPhiPtW);
 412     fHistList->Add(fh2TracksLeadingPhiPtW);
 413     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 414     fHistList->Add(fh2NRecJetsPtRan);
 415     fHistList->Add(fh2NConstPtRan);
 416     fHistList->Add(fh2NConstLeadingPtRan);
 417     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 418     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 419     }
 420
 421   // =========== Switch on Sumw2 for all histos ===========
 422   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 423     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 424     if (h1){
 425       h1->Sumw2();
 426       continue;
 427     }
 428     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 429     if(hn)hn->Sumw2();
 430   }
 431   TH1::AddDirectory(oldStatus);
 432 }
 433
 434 void AliAnalysisTaskJetCluster::Init()
 435 {
 436   //
 437   // Initialization
 438   //
 439
 440   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 441
 442 }
 443
 444 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 445 {
 446
 447   if(fUseGlobalSelection){
 448     // no selection by the service task, we continue
 449     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 450     PostData(1, fHistList);
 451     return;
 452   }
 453
 454   //
 455   // Execute analysis for current event
 456   //
 457   AliESDEvent *fESD = 0;
 458   if(fUseAODTrackInput){
 459     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 460     if(!fAOD){
 461       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 462       return;
 463     }
 464     // fethc the header
 465   }
 466   else{
 467     //  assume that the AOD is in the general output...
 468     fAOD  = AODEvent();
 469     if(!fAOD){
 470       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 471       return;
 472     }
 473     if(fDebug>0){
 474       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 475     }
 476   }
 477
 478   Bool_t selectEvent =  false;
 479   Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
 480
 481   if(fAOD){
 482     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 483     if(vtxAOD->GetNContributors()>0){
 484         Float_t zvtx = vtxAOD->GetZ();
 485         Float_t yvtx = vtxAOD->GetY();
 486         Float_t xvtx = vtxAOD->GetX();
 487         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 488         if(TMath::Abs(zvtx)<8.&&r2<1.){
 489           if(physicsSelection)selectEvent = true;
 490         }
 491     }
 492   }
 493   if(!selectEvent){
 494     PostData(1, fHistList);
 495     return;
 496   }
 497
 498   fh1Trials->Fill("#sum{ntrials}",1);
 499
 500
 501   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 502
 503   // ==== General variables needed
 504
 505
 506
 507   // we simply fetch the tracks/mc particles as a list of AliVParticles
 508
 509   TList recParticles;
 510   TList genParticles;
 511
 512   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 513   Float_t nCh = recParticles.GetEntries();
 514   fh1Nch->Fill(nCh);
 515   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 516   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 517   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 518
 519   // find the jets....
 520
 521   vector<fastjet::PseudoJet> inputParticlesRec;
 522   vector<fastjet::PseudoJet> inputParticlesRecRan;
 523   for(int i = 0; i < recParticles.GetEntries(); i++){
 524     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 525     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 526     // we talk total momentum here
 527     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 528     jInp.set_user_index(i);
 529     inputParticlesRec.push_back(jInp);
 530
 531     // the randomized input changes eta and phi, but keeps the p_T
 532     Double_t pT = vp->Pt();
 533     Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
 534     Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
 535
 536
 537     Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
 538     Double_t pZ = pT/TMath::Tan(theta);
 539
 540     Double_t pX = pT * TMath::Cos(phi);
 541     Double_t pY = pT * TMath::Sin(phi);
 542     Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 543
 544     fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 545     jInpRan.set_user_index(i);
 546     inputParticlesRecRan.push_back(jInpRan);
 547
 548   }
 549
 550   if(inputParticlesRec.size()==0){
 551     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 552     PostData(1, fHistList);
 553     return;
 554   }
 555
 556   // run fast jet
 557   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 558   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 559   fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
 560
 561   inclusiveJets = clustSeq.inclusive_jets();
 562   sortedJets    = sorted_by_pt(inclusiveJets);
 563
 564   fh1NJetsRec->Fill(sortedJets.size());
 565
 566  // loop over all jets an fill information, first on is the leading jet
 567
 568  Int_t nRecOver = inclusiveJets.size();
 569  Int_t nRec     = inclusiveJets.size();
 570  if(inclusiveJets.size()>0){
 571    AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 572     Float_t pt = leadingJet.Pt();
 573
 574     Int_t iCount = 0;
 575     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 576       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 577       while(pt<ptCut&&iCount<nRec){
 578         nRecOver--;
 579         iCount++;
 580         if(iCount<nRec){
 581           pt = sortedJets[iCount].perp();
 582         }
 583       }
 584       if(nRecOver<=0)break;
 585       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 586     }
 587     Float_t phi = leadingJet.Phi();
 588     if(phi<0)phi+=TMath::Pi()*2.;
 589     Float_t eta = leadingJet.Eta();
 590     pt = leadingJet.Pt();
 591
 592     // correlation of leading jet with tracks
 593     TIterator *recIter = recParticles.MakeIterator();
 594     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
 595
 596     recIter->Reset();
 597     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 598       Float_t tmpPt = tmpRecTrack->Pt();
 599       // correlation
 600       Float_t tmpPhi =  tmpRecTrack->Phi();
 601       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 602       Float_t dPhi = phi - tmpPhi;
 603       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 604       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 605       //      Float_t dEta = eta - tmpRecTrack->Eta();
 606       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 607       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 608    }
 609
 610
 611
 612     for(int j = 0; j < nRec;j++){
 613       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 614       Float_t tmpPt = tmpRec.Pt();
 615       fh1PtJetsRecIn->Fill(tmpPt);
 616       // Fill Spectra with constituents
 617       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 618       fh1NConstRec->Fill(constituents.size());
 619       fh2PtNch->Fill(nCh,tmpPt);
 620       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 621       fh2NConstPt->Fill(tmpPt,constituents.size());
 622       // loop over constiutents and fill spectrum
 623       for(unsigned int ic = 0; ic < constituents.size();ic++){
 624         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 625         fh1PtJetConstRec->Fill(part->Pt());
 626         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
 627       }
 628
 629       // correlation
 630       Float_t tmpPhi =  tmpRec.Phi();
 631       Float_t tmpEta =  tmpRec.Eta();
 632       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 633
 634       if(j==0){
 635         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 636         fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 637         fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 638         fh1NConstLeadingRec->Fill(constituents.size());
 639         fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 640         continue;
 641       }
 642       fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
 643       fh2JetEtaPt->Fill(tmpEta,tmpPt);
 644       Float_t dPhi = phi - tmpPhi;
 645       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 646       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 647       Float_t dEta = eta - tmpRec.Eta();
 648       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 649       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 650       fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 651     }
 652     delete recIter;
 653  }
 654
 655
 656  // fill track information
 657  Int_t nTrackOver = recParticles.GetSize();
 658   // do the same for tracks and jets
 659  if(nTrackOver>0){
 660    TIterator *recIter = recParticles.MakeIterator();
 661    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 662    Float_t pt = tmpRec->Pt();
 663    //    Printf("Leading track p_t %3.3E",pt);
 664    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 665      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 666      while(pt<ptCut&&tmpRec){
 667        nTrackOver--;
 668        tmpRec = (AliVParticle*)(recIter->Next());
 669        if(tmpRec){
 670          pt = tmpRec->Pt();
 671        }
 672      }
 673      if(nTrackOver<=0)break;
 674      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 675    }
 676
 677    recIter->Reset();
 678    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 679    Float_t phi = leading->Phi();
 680    if(phi<0)phi+=TMath::Pi()*2.;
 681    Float_t eta = leading->Eta();
 682    pt  = leading->Pt();
 683    while((tmpRec = (AliVParticle*)(recIter->Next()))){
 684      Float_t tmpPt = tmpRec->Pt();
 685      Float_t tmpEta = tmpRec->Eta();
 686      fh1PtTracksRecIn->Fill(tmpPt);
 687      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
 688      if(tmpRec==leading){
 689        fh1PtTracksLeadingRecIn->Fill(tmpPt);
 690        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
 691        continue;
 692      }
 693       // correlation
 694      Float_t tmpPhi =  tmpRec->Phi();
 695
 696      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 697      Float_t dPhi = phi - tmpPhi;
 698      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 699      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 700      Float_t dEta = eta - tmpRec->Eta();
 701      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 702      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 703      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 704    }
 705    delete recIter;
 706  }
 707
 708  // find the random jets
 709  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
 710  fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
 711
 712  inclusiveJetsRan = clustSeqRan.inclusive_jets();
 713  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
 714
 715  // fill the jet information from random track
 716
 717   fh1NJetsRecRan->Fill(sortedJetsRan.size());
 718
 719  // loop over all jets an fill information, first on is the leading jet
 720
 721  Int_t nRecOverRan = inclusiveJetsRan.size();
 722  Int_t nRecRan     = inclusiveJetsRan.size();
 723  if(inclusiveJetsRan.size()>0){
 724    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
 725    Float_t pt = leadingJet.Pt();
 726
 727    Int_t iCount = 0;
 728    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
 729      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
 730       while(pt<ptCut&&iCount<nRecRan){
 731         nRecOverRan--;
 732         iCount++;
 733         if(iCount<nRecRan){
 734           pt = sortedJetsRan[iCount].perp();
 735         }
 736       }
 737       if(nRecOverRan<=0)break;
 738       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
 739     }
 740     Float_t phi = leadingJet.Phi();
 741     if(phi<0)phi+=TMath::Pi()*2.;
 742     pt = leadingJet.Pt();
 743
 744     // correlation of leading jet with random tracks
 745
 746     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
 747       {
 748         Float_t tmpPt = inputParticlesRecRan[ip].perp();
 749         // correlation
 750         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
 751         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 752         Float_t dPhi = phi - tmpPhi;
 753         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 754         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 755         //      Float_t dEta = eta - tmpRecTrack->Eta();
 756         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
 757         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
 758       }
 759
 760
 761
 762     for(int j = 0; j < nRecRan;j++){
 763       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
 764       Float_t tmpPt = tmpRec.Pt();
 765       fh1PtJetsRecInRan->Fill(tmpPt);
 766       // Fill Spectra with constituents
 767       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
 768       fh1NConstRecRan->Fill(constituents.size());
 769       fh2NConstPtRan->Fill(tmpPt,constituents.size());
 770       fh2PtNchRan->Fill(nCh,tmpPt);
 771       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
 772       // correlation
 773       Float_t tmpPhi =  tmpRec.Phi();
 774       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 775
 776       if(j==0){
 777         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
 778         fh1NConstLeadingRecRan->Fill(constituents.size());
 779         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
 780         continue;
 781       }
 782     }
 783  }
 784
 785
 786  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 787  PostData(1, fHistList);
 788 }
 789
 790 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
 791 {
 792 // Terminate analysis
 793 //
 794     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 795 }
 796
 797
 798 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 799
 800   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 801
 802   Int_t iCount = 0;
 803   if(type==kTrackAOD){
 804     AliAODEvent *aod = 0;
 805     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 806     else aod = AODEvent();
 807     if(!aod){
 808       return iCount;
 809     }
 810     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 811       AliAODTrack *tr = aod->GetTrack(it);
 812       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 813       if(TMath::Abs(tr->Eta())>0.9)continue;
 814       //      if(tr->Pt()<0.3)continue;
 815       list->Add(tr);
 816       iCount++;
 817     }
 818   }
 819   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 820     AliMCEvent* mcEvent = MCEvent();
 821     if(!mcEvent)return iCount;
 822     // we want to have alivpartilces so use get track
 823     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 824       if(!mcEvent->IsPhysicalPrimary(it))continue;
 825       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 826       if(type == kTrackKineAll){
 827         list->Add(part);
 828         iCount++;
 829       }
 830       else if(type == kTrackKineCharged){
 831         if(part->Particle()->GetPDG()->Charge()==0)continue;
 832         list->Add(part);
 833         iCount++;
 834       }
 835     }
 836   }
 837   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 838     AliAODEvent *aod = 0;
 839     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 840     else aod = AODEvent();
 841     if(!aod)return iCount;
 842     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 843     if(!tca)return iCount;
 844     for(int it = 0;it < tca->GetEntriesFast();++it){
 845       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 846       if(!part->IsPhysicalPrimary())continue;
 847       if(type == kTrackAODMCAll){
 848         list->Add(part);
 849         iCount++;
 850       }
 851       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 852         if(part->Charge()==0)continue;
 853         if(kTrackAODMCCharged){
 854           list->Add(part);
 855         }
 856         else {
 857           if(TMath::Abs(part->Eta())>0.9)continue;
 858           list->Add(part);
 859         }
 860         iCount++;
 861       }
 862     }
 863   }// AODMCparticle
 864   list->Sort();
 865   return iCount;
 866
 867 }
 868
 869 /*
 870 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
 871   for(int i = 0; i < particles.GetEntries(); i++){
 872     AliVParticle *vp = (AliVParticle*)particles.At(i);
 873     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 874     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 875     jInp.set_user_index(i);
 876     inputParticles.push_back(jInp);
 877   }
 878
 879   return 0;
 880
 881 }
 882 */