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   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 417
 418 }
 419
 420 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 421 {
 422
 423   if(fUseGlobalSelection){
 424     // no selection by the service task, we continue
 425     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 426     PostData(1, fHistList);
 427     return;
 428   }
 429
 430   //
 431   // Execute analysis for current event
 432   //
 433   AliESDEvent *fESD = 0;
 434   if(fUseAODTrackInput){
 435     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 436     if(!fAOD){
 437       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 438       return;
 439     }
 440     // fethc the header
 441   }
 442   else{
 443     //  assume that the AOD is in the general output...
 444     fAOD  = AODEvent();
 445     if(!fAOD){
 446       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 447       return;
 448     }
 449     if(fDebug>0){
 450       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 451     }
 452   }
 453
 454   Bool_t selectEvent =  false;
 455   Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
 456
 457   if(fAOD){
 458     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 459     if(vtxAOD->GetNContributors()>0){
 460         Float_t zvtx = vtxAOD->GetZ();
 461         Float_t yvtx = vtxAOD->GetY();
 462         Float_t xvtx = vtxAOD->GetX();
 463         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 464         if(TMath::Abs(zvtx)<8.&&r2<1.){
 465           if(physicsSelection)selectEvent = true;
 466         }
 467     }
 468   }
 469   if(!selectEvent){
 470     PostData(1, fHistList);
 471     return;
 472   }
 473
 474   fh1Trials->Fill("#sum{ntrials}",1);
 475
 476
 477   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 478
 479   // ==== General variables needed
 480
 481
 482
 483   // we simply fetch the tracks/mc particles as a list of AliVParticles
 484
 485   TList recParticles;
 486   TList genParticles;
 487
 488   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 489   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 490   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 491   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 492
 493   // find the jets....
 494
 495   vector<fastjet::PseudoJet> inputParticlesRec;
 496   vector<fastjet::PseudoJet> inputParticlesRecRan;
 497   for(int i = 0; i < recParticles.GetEntries(); i++){
 498     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 499     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 500     // we talk total momentum here
 501     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 502     jInp.set_user_index(i);
 503     inputParticlesRec.push_back(jInp);
 504
 505     // the randomized input changes eta and phi, but keeps the p_T
 506     Double_t pT = vp->Pt();
 507     Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
 508     Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
 509
 510
 511     Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
 512     Double_t pZ = pT/TMath::Tan(theta);
 513
 514     Double_t pX = pT * TMath::Cos(phi);
 515     Double_t pY = pT * TMath::Sin(phi);
 516     Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 517
 518     fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 519     jInpRan.set_user_index(i);
 520     inputParticlesRecRan.push_back(jInpRan);
 521
 522   }
 523
 524   if(inputParticlesRec.size()==0){
 525     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 526     PostData(1, fHistList);
 527     return;
 528   }
 529
 530   // run fast jet
 531   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 532   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 533   fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
 534
 535   inclusiveJets = clustSeq.inclusive_jets();
 536   sortedJets    = sorted_by_pt(inclusiveJets);
 537
 538   if(fDebug)Printf("%s:%d Found %d jets",(char*)__FILE__,__LINE__, inclusiveJets.size());
 539
 540   fh1NJetsRec->Fill(sortedJets.size());
 541
 542  // loop over all jets an fill information, first on is the leading jet
 543
 544  Int_t nRecOver = inclusiveJets.size();
 545  Int_t nRec     = inclusiveJets.size();
 546  if(inclusiveJets.size()>0){
 547    AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 548     Float_t pt = leadingJet.Pt();
 549
 550     Int_t iCount = 0;
 551     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 552       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 553       while(pt<ptCut&&iCount<nRec){
 554         nRecOver--;
 555         iCount++;
 556         if(iCount<nRec){
 557           pt = sortedJets[iCount].perp();
 558         }
 559       }
 560       if(nRecOver<=0)break;
 561       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 562     }
 563     Float_t phi = leadingJet.Phi();
 564     if(phi<0)phi+=TMath::Pi()*2.;
 565     Float_t eta = leadingJet.Eta();
 566     pt = leadingJet.Pt();
 567
 568     // correlation of leading jet with tracks
 569     TIterator *recIter = recParticles.MakeIterator();
 570     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
 571
 572     recIter->Reset();
 573     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 574       Float_t tmpPt = tmpRecTrack->Pt();
 575       // correlation
 576       Float_t tmpPhi =  tmpRecTrack->Phi();
 577       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 578       Float_t dPhi = phi - tmpPhi;
 579       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 580       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 581       //      Float_t dEta = eta - tmpRecTrack->Eta();
 582       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 583       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 584    }
 585
 586
 587
 588     for(int j = 0; j < nRec;j++){
 589       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 590       Float_t tmpPt = tmpRec.Pt();
 591       fh1PtJetsRecIn->Fill(tmpPt);
 592       // Fill Spectra with constituents
 593       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 594       fh1NConstRec->Fill(constituents.size());
 595       fh2NConstPt->Fill(tmpPt,constituents.size());
 596       // loop over constiutents and fill spectrum
 597       for(unsigned int ic = 0; ic < constituents.size();ic++){
 598         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 599         fh1PtJetConstRec->Fill(part->Pt());
 600         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
 601       }
 602
 603       // correlation
 604       Float_t tmpPhi =  tmpRec.Phi();
 605       Float_t tmpEta =  tmpRec.Eta();
 606       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 607
 608       if(j==0){
 609         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 610         fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 611         fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 612         fh1NConstLeadingRec->Fill(constituents.size());
 613         fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 614         continue;
 615       }
 616       fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
 617       fh2JetEtaPt->Fill(tmpEta,tmpPt);
 618       Float_t dPhi = phi - tmpPhi;
 619       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 620       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 621       Float_t dEta = eta - tmpRec.Eta();
 622       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 623       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 624       fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 625     }
 626     delete recIter;
 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 */