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