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   if(fDebug)Printf("%s:%d Found %d jets",(char*)__FILE__,__LINE__, inclusiveJets.size());
 565
 566   fh1NJetsRec->Fill(sortedJets.size());
 567
 568  // loop over all jets an fill information, first on is the leading jet
 569
 570  Int_t nRecOver = inclusiveJets.size();
 571  Int_t nRec     = inclusiveJets.size();
 572  if(inclusiveJets.size()>0){
 573    AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 574     Float_t pt = leadingJet.Pt();
 575
 576     Int_t iCount = 0;
 577     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 578       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 579       while(pt<ptCut&&iCount<nRec){
 580         nRecOver--;
 581         iCount++;
 582         if(iCount<nRec){
 583           pt = sortedJets[iCount].perp();
 584         }
 585       }
 586       if(nRecOver<=0)break;
 587       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 588     }
 589     Float_t phi = leadingJet.Phi();
 590     if(phi<0)phi+=TMath::Pi()*2.;
 591     Float_t eta = leadingJet.Eta();
 592     pt = leadingJet.Pt();
 593
 594     // correlation of leading jet with tracks
 595     TIterator *recIter = recParticles.MakeIterator();
 596     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
 597
 598     recIter->Reset();
 599     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 600       Float_t tmpPt = tmpRecTrack->Pt();
 601       // correlation
 602       Float_t tmpPhi =  tmpRecTrack->Phi();
 603       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 604       Float_t dPhi = phi - tmpPhi;
 605       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 606       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 607       //      Float_t dEta = eta - tmpRecTrack->Eta();
 608       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 609       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 610    }
 611
 612
 613
 614     for(int j = 0; j < nRec;j++){
 615       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 616       Float_t tmpPt = tmpRec.Pt();
 617       fh1PtJetsRecIn->Fill(tmpPt);
 618       // Fill Spectra with constituents
 619       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 620       fh1NConstRec->Fill(constituents.size());
 621       fh2PtNch->Fill(nCh,tmpPt);
 622       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 623       fh2NConstPt->Fill(tmpPt,constituents.size());
 624       // loop over constiutents and fill spectrum
 625       for(unsigned int ic = 0; ic < constituents.size();ic++){
 626         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 627         fh1PtJetConstRec->Fill(part->Pt());
 628         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
 629       }
 630
 631       // correlation
 632       Float_t tmpPhi =  tmpRec.Phi();
 633       Float_t tmpEta =  tmpRec.Eta();
 634       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 635
 636       if(j==0){
 637         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 638         fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 639         fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 640         fh1NConstLeadingRec->Fill(constituents.size());
 641         fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 642         continue;
 643       }
 644       fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
 645       fh2JetEtaPt->Fill(tmpEta,tmpPt);
 646       Float_t dPhi = phi - tmpPhi;
 647       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 648       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 649       Float_t dEta = eta - tmpRec.Eta();
 650       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 651       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 652       fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 653     }
 654     delete recIter;
 655  }
 656
 657
 658  // fill track information
 659  Int_t nTrackOver = recParticles.GetSize();
 660   // do the same for tracks and jets
 661  if(nTrackOver>0){
 662    TIterator *recIter = recParticles.MakeIterator();
 663    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 664    Float_t pt = tmpRec->Pt();
 665    //    Printf("Leading track p_t %3.3E",pt);
 666    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 667      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 668      while(pt<ptCut&&tmpRec){
 669        nTrackOver--;
 670        tmpRec = (AliVParticle*)(recIter->Next());
 671        if(tmpRec){
 672          pt = tmpRec->Pt();
 673        }
 674      }
 675      if(nTrackOver<=0)break;
 676      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 677    }
 678
 679    recIter->Reset();
 680    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 681    Float_t phi = leading->Phi();
 682    if(phi<0)phi+=TMath::Pi()*2.;
 683    Float_t eta = leading->Eta();
 684    pt  = leading->Pt();
 685    while((tmpRec = (AliVParticle*)(recIter->Next()))){
 686      Float_t tmpPt = tmpRec->Pt();
 687      Float_t tmpEta = tmpRec->Eta();
 688      fh1PtTracksRecIn->Fill(tmpPt);
 689      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
 690      if(tmpRec==leading){
 691        fh1PtTracksLeadingRecIn->Fill(tmpPt);
 692        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
 693        continue;
 694      }
 695       // correlation
 696      Float_t tmpPhi =  tmpRec->Phi();
 697
 698      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 699      Float_t dPhi = phi - tmpPhi;
 700      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 701      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 702      Float_t dEta = eta - tmpRec->Eta();
 703      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 704      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 705      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 706    }
 707    delete recIter;
 708  }
 709
 710  // find the random jets
 711  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
 712  fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
 713
 714  inclusiveJetsRan = clustSeqRan.inclusive_jets();
 715  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
 716
 717  // fill the jet information from random track
 718
 719   fh1NJetsRecRan->Fill(sortedJetsRan.size());
 720
 721  // loop over all jets an fill information, first on is the leading jet
 722
 723  Int_t nRecOverRan = inclusiveJetsRan.size();
 724  Int_t nRecRan     = inclusiveJetsRan.size();
 725  if(inclusiveJetsRan.size()>0){
 726    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
 727    Float_t pt = leadingJet.Pt();
 728
 729    Int_t iCount = 0;
 730    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
 731      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
 732       while(pt<ptCut&&iCount<nRecRan){
 733         nRecOverRan--;
 734         iCount++;
 735         if(iCount<nRecRan){
 736           pt = sortedJetsRan[iCount].perp();
 737         }
 738       }
 739       if(nRecOverRan<=0)break;
 740       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
 741     }
 742     Float_t phi = leadingJet.Phi();
 743     if(phi<0)phi+=TMath::Pi()*2.;
 744     pt = leadingJet.Pt();
 745
 746     // correlation of leading jet with random tracks
 747
 748     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
 749       {
 750         Float_t tmpPt = inputParticlesRecRan[ip].perp();
 751         // correlation
 752         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
 753         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 754         Float_t dPhi = phi - tmpPhi;
 755         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 756         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 757         //      Float_t dEta = eta - tmpRecTrack->Eta();
 758         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
 759         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
 760       }
 761
 762
 763
 764     for(int j = 0; j < nRecRan;j++){
 765       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
 766       Float_t tmpPt = tmpRec.Pt();
 767       fh1PtJetsRecInRan->Fill(tmpPt);
 768       // Fill Spectra with constituents
 769       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
 770       fh1NConstRecRan->Fill(constituents.size());
 771       fh2NConstPtRan->Fill(tmpPt,constituents.size());
 772       fh2PtNchRan->Fill(nCh,tmpPt);
 773       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
 774       // correlation
 775       Float_t tmpPhi =  tmpRec.Phi();
 776       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 777
 778       if(j==0){
 779         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
 780         fh1NConstLeadingRecRan->Fill(constituents.size());
 781         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
 782         continue;
 783       }
 784     }
 785  }
 786
 787
 788  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 789  PostData(1, fHistList);
 790 }
 791
 792 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
 793 {
 794 // Terminate analysis
 795 //
 796     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 797 }
 798
 799
 800 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 801
 802   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 803
 804   Int_t iCount = 0;
 805   if(type==kTrackAOD){
 806     AliAODEvent *aod = 0;
 807     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 808     else aod = AODEvent();
 809     if(!aod){
 810       return iCount;
 811     }
 812     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 813       AliAODTrack *tr = aod->GetTrack(it);
 814       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 815       if(TMath::Abs(tr->Eta())>0.9)continue;
 816       //      if(tr->Pt()<0.3)continue;
 817       list->Add(tr);
 818       iCount++;
 819     }
 820   }
 821   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 822     AliMCEvent* mcEvent = MCEvent();
 823     if(!mcEvent)return iCount;
 824     // we want to have alivpartilces so use get track
 825     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 826       if(!mcEvent->IsPhysicalPrimary(it))continue;
 827       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 828       if(type == kTrackKineAll){
 829         list->Add(part);
 830         iCount++;
 831       }
 832       else if(type == kTrackKineCharged){
 833         if(part->Particle()->GetPDG()->Charge()==0)continue;
 834         list->Add(part);
 835         iCount++;
 836       }
 837     }
 838   }
 839   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 840     AliAODEvent *aod = 0;
 841     if(fUseAODMCInput)dynamic_cast<AliAODEvent*>(InputEvent());
 842     else aod = AODEvent();
 843     if(!aod)return iCount;
 844     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 845     if(!tca)return iCount;
 846     for(int it = 0;it < tca->GetEntriesFast();++it){
 847       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 848       if(!part->IsPhysicalPrimary())continue;
 849       if(type == kTrackAODMCAll){
 850         list->Add(part);
 851         iCount++;
 852       }
 853       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 854         if(part->Charge()==0)continue;
 855         if(kTrackAODMCCharged){
 856           list->Add(part);
 857         }
 858         else {
 859           if(TMath::Abs(part->Eta())>0.9)continue;
 860           list->Add(part);
 861         }
 862         iCount++;
 863       }
 864     }
 865   }// AODMCparticle
 866   list->Sort();
 867   return iCount;
 868
 869 }
 870
 871 /*
 872 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
 873   for(int i = 0; i < particles.GetEntries(); i++){
 874     AliVParticle *vp = (AliVParticle*)particles.At(i);
 875     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 876     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 877     jInp.set_user_index(i);
 878     inputParticles.push_back(jInp);
 879   }
 880
 881   return 0;
 882
 883 }
 884 */