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     delete recIter;
 628  }
 629
 630
 631  // fill track information
 632  Int_t nTrackOver = recParticles.GetSize();
 633   // do the same for tracks and jets
 634  if(nTrackOver>0){
 635    TIterator *recIter = recParticles.MakeIterator();
 636    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 637    Float_t pt = tmpRec->Pt();
 638    //    Printf("Leading track p_t %3.3E",pt);
 639    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 640      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 641      while(pt<ptCut&&tmpRec){
 642        nTrackOver--;
 643        tmpRec = (AliVParticle*)(recIter->Next());
 644        if(tmpRec){
 645          pt = tmpRec->Pt();
 646        }
 647      }
 648      if(nTrackOver<=0)break;
 649      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 650    }
 651
 652    recIter->Reset();
 653    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 654    Float_t phi = leading->Phi();
 655    if(phi<0)phi+=TMath::Pi()*2.;
 656    Float_t eta = leading->Eta();
 657    pt  = leading->Pt();
 658    while((tmpRec = (AliVParticle*)(recIter->Next()))){
 659      Float_t tmpPt = tmpRec->Pt();
 660      Float_t tmpEta = tmpRec->Eta();
 661      fh1PtTracksRecIn->Fill(tmpPt);
 662      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
 663      if(tmpRec==leading){
 664        fh1PtTracksLeadingRecIn->Fill(tmpPt);
 665        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
 666        continue;
 667      }
 668       // correlation
 669      Float_t tmpPhi =  tmpRec->Phi();
 670
 671      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 672      Float_t dPhi = phi - tmpPhi;
 673      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 674      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 675      Float_t dEta = eta - tmpRec->Eta();
 676      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 677      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 678      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 679    }
 680    delete recIter;
 681  }
 682
 683  // find the random jets
 684  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
 685  fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
 686
 687  inclusiveJetsRan = clustSeqRan.inclusive_jets();
 688  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
 689
 690  // fill the jet information from random track
 691
 692   fh1NJetsRecRan->Fill(sortedJetsRan.size());
 693
 694  // loop over all jets an fill information, first on is the leading jet
 695
 696  Int_t nRecOverRan = inclusiveJetsRan.size();
 697  Int_t nRecRan     = inclusiveJetsRan.size();
 698  if(inclusiveJetsRan.size()>0){
 699    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
 700    Float_t pt = leadingJet.Pt();
 701
 702    Int_t iCount = 0;
 703    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
 704      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
 705       while(pt<ptCut&&iCount<nRecRan){
 706         nRecOverRan--;
 707         iCount++;
 708         if(iCount<nRecRan){
 709           pt = sortedJetsRan[iCount].perp();
 710         }
 711       }
 712       if(nRecOverRan<=0)break;
 713       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
 714     }
 715     Float_t phi = leadingJet.Phi();
 716     if(phi<0)phi+=TMath::Pi()*2.;
 717     pt = leadingJet.Pt();
 718
 719     // correlation of leading jet with random tracks
 720
 721     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
 722       {
 723         Float_t tmpPt = inputParticlesRecRan[ip].perp();
 724         // correlation
 725         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
 726         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 727         Float_t dPhi = phi - tmpPhi;
 728         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 729         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 730         //      Float_t dEta = eta - tmpRecTrack->Eta();
 731         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
 732         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
 733       }
 734
 735
 736
 737     for(int j = 0; j < nRecRan;j++){
 738       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
 739       Float_t tmpPt = tmpRec.Pt();
 740       fh1PtJetsRecInRan->Fill(tmpPt);
 741       // Fill Spectra with constituents
 742       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
 743       fh1NConstRecRan->Fill(constituents.size());
 744       fh2NConstPtRan->Fill(tmpPt,constituents.size());
 745
 746       // correlation
 747       Float_t tmpPhi =  tmpRec.Phi();
 748       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 749
 750       if(j==0){
 751         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
 752         fh1NConstLeadingRecRan->Fill(constituents.size());
 753         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
 754         continue;
 755       }
 756     }
 757  }
 758
 759
 760  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 761  PostData(1, fHistList);
 762 }
 763
 764 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
 765 {
 766 // Terminate analysis
 767 //
 768     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 769 }
 770
 771
 772 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 773
 774   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 775
 776   Int_t iCount = 0;
 777   if(type==kTrackAOD){
 778     AliAODEvent *aod = 0;
 779     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 780     else aod = AODEvent();
 781     if(!aod){
 782       return iCount;
 783     }
 784     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 785       AliAODTrack *tr = aod->GetTrack(it);
 786       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 787       if(TMath::Abs(tr->Eta())>0.9)continue;
 788       //      if(tr->Pt()<0.3)continue;
 789       list->Add(tr);
 790       iCount++;
 791     }
 792   }
 793   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 794     AliMCEvent* mcEvent = MCEvent();
 795     if(!mcEvent)return iCount;
 796     // we want to have alivpartilces so use get track
 797     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 798       if(!mcEvent->IsPhysicalPrimary(it))continue;
 799       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 800       if(type == kTrackKineAll){
 801         list->Add(part);
 802         iCount++;
 803       }
 804       else if(type == kTrackKineCharged){
 805         if(part->Particle()->GetPDG()->Charge()==0)continue;
 806         list->Add(part);
 807         iCount++;
 808       }
 809     }
 810   }
 811   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 812     AliAODEvent *aod = 0;
 813     if(fUseAODMCInput)dynamic_cast<AliAODEvent*>(InputEvent());
 814     else aod = AODEvent();
 815     if(!aod)return iCount;
 816     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 817     if(!tca)return iCount;
 818     for(int it = 0;it < tca->GetEntriesFast();++it){
 819       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 820       if(!part->IsPhysicalPrimary())continue;
 821       if(type == kTrackAODMCAll){
 822         list->Add(part);
 823         iCount++;
 824       }
 825       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 826         if(part->Charge()==0)continue;
 827         if(kTrackAODMCCharged){
 828           list->Add(part);
 829         }
 830         else {
 831           if(TMath::Abs(part->Eta())>0.9)continue;
 832           list->Add(part);
 833         }
 834         iCount++;
 835       }
 836     }
 837   }// AODMCparticle
 838   list->Sort();
 839   return iCount;
 840
 841 }
 842
 843 /*
 844 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
 845   for(int i = 0; i < particles.GetEntries(); i++){
 846     AliVParticle *vp = (AliVParticle*)particles.At(i);
 847     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 848     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 849     jInp.set_user_index(i);
 850     inputParticles.push_back(jInp);
 851   }
 852
 853   return 0;
 854
 855 }
 856 */