PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

   1 // **************************************
   2 //  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 <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TClonesArray.h>
  37 #include <TRefArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetSpectrum2.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 "AliAODJetEventBackground.h"
  53 #include "AliAODMCParticle.h"
  54 #include "AliMCEventHandler.h"
  55 #include "AliMCEvent.h"
  56 #include "AliStack.h"
  57 #include "AliGenPythiaEventHeader.h"
  58 #include "AliJetKineReaderHeader.h"
  59 #include "AliGenCocktailEventHeader.h"
  60 #include "AliInputEventHandler.h"
  61
  62
  63 #include "AliAnalysisHelperJetTasks.h"
  64
  65 ClassImp(AliAnalysisTaskJetSpectrum2)
  66
  67 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
  68   AliAnalysisTaskSE(),
  69   fJetHeaderRec(0x0),
  70   fJetHeaderGen(0x0),
  71   fAODIn(0x0),
  72   fAODOut(0x0),
  73   fAODExtension(0x0),
  74   fhnCorrelation(0x0),
  75   fhnCorrelationPhiZRec(0x0),
  76   f1PtScale(0x0),
  77   fBranchRec("jets"),
  78   fBranchGen(""),
  79   fBranchBkgRec(""),
  80   fBranchBkgGen(""),
  81   fNonStdFile(""),
  82   fUseAODJetInput(kFALSE),
  83   fUseAODTrackInput(kFALSE),
  84   fUseAODMCInput(kFALSE),
  85   fUseGlobalSelection(kFALSE),
  86   fUseExternalWeightOnly(kFALSE),
  87   fLimitGenJetEta(kFALSE),
  88   fNMatchJets(5),
  89   fFilterMask(0),
  90   fEventSelectionMask(0),
  91   fAnalysisType(0),
  92   fTrackTypeRec(kTrackUndef),
  93   fTrackTypeGen(kTrackUndef),
  94   fEventClass(0),
  95   fAvgTrials(1),
  96   fExternalWeight(1),
  97   fJetRecEtaWindow(0.5),
  98   fTrackRecEtaWindow(0.5),
  99   fMinJetPt(0),
 100   fMinTrackPt(0.15),
 101   fDeltaPhiWindow(90./180.*TMath::Pi()),
 102   fMultRec(0),
 103   fMultGen(0),
 104   fh1Xsec(0x0),
 105   fh1Trials(0x0),
 106   fh1PtHard(0x0),
 107   fh1PtHardNoW(0x0),
 108   fh1PtHardTrials(0x0),
 109   fh1ZVtx(0x0),
 110   fh1TmpRho(0x0),
 111   fh2MultRec(0x0),
 112   fh2MultGen(0x0),
 113   fh2PtFGen(0x0),
 114   fh2RelPtFGen(0x0),
 115   fHistList(0x0)
 116 {
 117   for(int i = 0;i < kMaxStep*2;++i){
 118     fhnJetContainer[i] = 0;
 119   }
 120
 121   for(int ij = 0;ij <kJetTypes;++ij){
 122     fh1NJets[ij] = 0;
 123     fh1SumPtTrack[ij] = 0;
 124     fh1PtJetsIn[ij] = 0;
 125     fh1PtTracksIn[ij] = 0;
 126     fh1PtTracksLeadingIn[ij] = 0;
 127     fh2MultJetPt[ij] = 0;
 128     fh2NJetsPt[ij]  = 0;
 129     fh2NTracksPt[ij]  = 0;
 130     fh2LeadingTrackPtTrackPhi[ij] = 0;
 131     for(int i = 0;i <= kMaxJets;++i){
 132       fh2PhiPt[ij][i] = 0;
 133       fh2EtaPt[ij][i] = 0;
 134       fh2AreaPt[ij][i] = 0;
 135       fh2EtaArea[ij][i] = 0;
 136       fh2PhiEta[ij][i] = 0;
 137       fh2LTrackPtJetPt[ij][i] = 0;
 138       fh1PtIn[ij][i] = 0;
 139       fh2RhoPt[ij][i] = 0;
 140       fh2PsiPt[ij][i] = 0;
 141     }
 142
 143     fh1DijetMinv[ij] = 0;
 144     fh2DijetDeltaPhiPt[ij] = 0;
 145     fh2DijetAsymPt[ij] = 0;
 146     fh2DijetPt2vsPt1[ij] = 0;
 147     fh2DijetDifvsSum[ij] = 0;
 148
 149   }
 150 }
 151
 152 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 153   AliAnalysisTaskSE(name),
 154   fJetHeaderRec(0x0),
 155   fJetHeaderGen(0x0),
 156   fAODIn(0x0),
 157   fAODOut(0x0),
 158   fAODExtension(0x0),
 159   fhnCorrelation(0x0),
 160   fhnCorrelationPhiZRec(0x0),
 161   f1PtScale(0x0),
 162   fBranchRec("jets"),
 163   fBranchGen(""),
 164   fBranchBkgRec(""),
 165   fBranchBkgGen(""),
 166   fNonStdFile(""),
 167   fUseAODJetInput(kFALSE),
 168   fUseAODTrackInput(kFALSE),
 169   fUseAODMCInput(kFALSE),
 170   fUseGlobalSelection(kFALSE),
 171   fUseExternalWeightOnly(kFALSE),
 172   fLimitGenJetEta(kFALSE),
 173   fNMatchJets(5),
 174   fFilterMask(0),
 175   fEventSelectionMask(0),
 176   fAnalysisType(0),
 177   fTrackTypeRec(kTrackUndef),
 178   fTrackTypeGen(kTrackUndef),
 179   fEventClass(0),
 180   fAvgTrials(1),
 181   fExternalWeight(1),
 182   fJetRecEtaWindow(0.5),
 183   fTrackRecEtaWindow(0.5),
 184   fMinJetPt(0),
 185   fMinTrackPt(0.15),
 186   fDeltaPhiWindow(90./180.*TMath::Pi()),
 187   fMultRec(0),
 188   fMultGen(0),
 189   fh1Xsec(0x0),
 190   fh1Trials(0x0),
 191   fh1PtHard(0x0),
 192   fh1PtHardNoW(0x0),
 193   fh1PtHardTrials(0x0),
 194   fh1ZVtx(0x0),
 195   fh1TmpRho(0x0),
 196   fh2MultRec(0x0),
 197   fh2MultGen(0x0),
 198   fh2PtFGen(0x0),
 199   fh2RelPtFGen(0x0),
 200   fHistList(0x0)
 201 {
 202
 203   for(int i = 0;i < kMaxStep*2;++i){
 204     fhnJetContainer[i] = 0;
 205   }
 206
 207   for(int ij = 0;ij <kJetTypes;++ij){
 208     fh1NJets[ij] = 0;
 209     fh1SumPtTrack[ij] = 0;
 210     fh1PtJetsIn[ij] = 0;
 211     fh1PtTracksIn[ij] = 0;
 212     fh1PtTracksLeadingIn[ij] = 0;
 213     fh2MultJetPt[ij] = 0;
 214     fh2NJetsPt[ij]  = 0;
 215     fh2NTracksPt[ij]  = 0;
 216     fh2LeadingTrackPtTrackPhi[ij] = 0;
 217     for(int i = 0;i <= kMaxJets;++i){
 218       fh2PhiPt[ij][i] = 0;
 219       fh2EtaPt[ij][i] = 0;
 220       fh2AreaPt[ij][i] = 0;
 221       fh2EtaArea[ij][i] = 0;
 222       fh2PhiEta[ij][i] = 0;
 223       fh2LTrackPtJetPt[ij][i] = 0;
 224       fh1PtIn[ij][i] = 0;
 225       fh2RhoPt[ij][i] = 0;
 226       fh2PsiPt[ij][i] = 0;
 227     }
 228
 229     fh1DijetMinv[ij] = 0;
 230     fh2DijetDeltaPhiPt[ij] = 0;
 231     fh2DijetAsymPt[ij] = 0;
 232     fh2DijetPt2vsPt1[ij] = 0;
 233     fh2DijetDifvsSum[ij] = 0;
 234   }
 235
 236   DefineOutput(1, TList::Class());
 237 }
 238
 239
 240
 241 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 242 {
 243
 244
 245
 246   //
 247   // Implemented Notify() to read the cross sections
 248   // and number of trials from pyxsec.root
 249   //
 250
 251   // Fetch the aod also from the input in,
 252   // have todo it in notify
 253
 254
 255   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 256   //  assume that the AOD is in the general output...
 257   fAODOut  = AODEvent();
 258
 259   if(fNonStdFile.Length()!=0){
 260     // case that we have an AOD extension we need can fetch the jets from the extended output
 261     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 262     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 263     if(!fAODExtension){
 264       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 265     }
 266   }
 267
 268
 269   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 270   Float_t xsection = 0;
 271   Float_t ftrials  = 1;
 272
 273   fAvgTrials = 1;
 274   if(tree){
 275     TFile *curfile = tree->GetCurrentFile();
 276     if (!curfile) {
 277       Error("Notify","No current file");
 278       return kFALSE;
 279     }
 280     if(!fh1Xsec||!fh1Trials){
 281       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 282       return kFALSE;
 283     }
 284     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 285     fh1Xsec->Fill("<#sigma>",xsection);
 286     // construct a poor man average trials
 287     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 288     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 289   }
 290
 291   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 292
 293   return kTRUE;
 294 }
 295
 296 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 297 {
 298
 299
 300   // Connect the AOD
 301
 302   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 303   OpenFile(1);
 304   if(!fHistList)fHistList = new TList();
 305   PostData(1, fHistList); // post data in any case once
 306
 307   fHistList->SetOwner(kTRUE);
 308   Bool_t oldStatus = TH1::AddDirectoryStatus();
 309   TH1::AddDirectory(kFALSE);
 310
 311   MakeJetContainer();
 312   fHistList->Add(fhnCorrelation);
 313   fHistList->Add(fhnCorrelationPhiZRec);
 314   for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 315
 316   //
 317   //  Histogram
 318
 319   const Int_t nBinPt = 320;
 320   Double_t binLimitsPt[nBinPt+1];
 321   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 322     if(iPt == 0){
 323       binLimitsPt[iPt] = 0.0;
 324     }
 325     else {// 1.0
 326       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 327     }
 328   }
 329   const Int_t nBinPhi = 90;
 330   Double_t binLimitsPhi[nBinPhi+1];
 331   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 332     if(iPhi==0){
 333       binLimitsPhi[iPhi] = 0;
 334     }
 335     else{
 336       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 337     }
 338   }
 339
 340   const Int_t nBinEta = 40;
 341   Double_t binLimitsEta[nBinEta+1];
 342   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 343     if(iEta==0){
 344       binLimitsEta[iEta] = -2.0;
 345     }
 346     else{
 347       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 348     }
 349   }
 350
 351
 352   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 353   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 354    fHistList->Add(fh1Xsec);
 355   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 356   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 357   fHistList->Add(fh1Trials);
 358   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 359   fHistList->Add(fh1PtHard);
 360   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 361   fHistList->Add(fh1PtHardNoW);
 362   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 363   fHistList->Add(fh1PtHardTrials);
 364
 365   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 366   fHistList->Add(fh1ZVtx);
 367   fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
 368   fHistList->Add(fh2MultRec);
 369   fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
 370   fHistList->Add(fh2MultGen);
 371
 372   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 373   fHistList->Add(fh2PtFGen);
 374
 375   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
 376   fHistList->Add(fh2RelPtFGen);
 377
 378     for(int ij = 0;ij <kJetTypes;++ij){
 379       TString cAdd = "";
 380       TString cJetBranch = "";
 381       if(ij==kJetRec){
 382         cAdd = "Rec";
 383         cJetBranch = fBranchRec.Data();
 384       }
 385       else if (ij==kJetGen){
 386         cAdd = "Gen";
 387         cJetBranch = fBranchGen.Data();
 388       }
 389       else if (ij==kJetRecFull){
 390         cAdd = "RecFull";
 391         cJetBranch = fBranchRec.Data();
 392       }
 393       else if (ij==kJetGenFull){
 394         cAdd = "GenFull";
 395         cJetBranch = fBranchGen.Data();
 396       }
 397
 398       fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
 399       fHistList->Add(fh1NJets[ij]);
 400
 401       fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 402       fHistList->Add(fh1PtJetsIn[ij]);
 403
 404       fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 405       fHistList->Add(fh1PtTracksIn[ij]);
 406
 407       fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 408       fHistList->Add(fh1PtTracksLeadingIn[ij]);
 409
 410       fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 411       fHistList->Add(fh1SumPtTrack[ij]);
 412
 413       fh2MultJetPt[ij]  = new TH2F(Form("fh2MultJetPt%s",cAdd.Data()),Form("%s jets p_T;# tracks;;p_{T} (GeV/c)",cAdd.Data()),400,0,4000,nBinPt,binLimitsPt);
 414       fHistList->Add(fh2MultJetPt[ij]);
 415
 416       fh2NJetsPt[ij]  = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
 417       fHistList->Add(fh2NJetsPt[ij]);
 418
 419       fh2NTracksPt[ij]  = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,4000);
 420       fHistList->Add(fh2NTracksPt[ij]);
 421
 422       fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
 423                                                nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
 424       fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
 425
 426       for(int i = 0;i <= kMaxJets;++i){
 427         fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 428         fHistList->Add(fh1PtIn[ij][i]);
 429
 430         fh2RhoPt[ij][i] = new TH2F(Form("fh2RhoPt%s_j%d",cAdd.Data(),i),Form("jet shape rho for %s jets;r;p_{T};",cAdd.Data()),
 431                                    40,0.,2.,nBinPt,binLimitsPt);
 432         fHistList->Add(fh2RhoPt[ij][i]);
 433         if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 434         fh2PsiPt[ij][i] = new TH2F(Form("fh2PsiPt%s_j%d",cAdd.Data(),i),Form("jet shape #psi for %s jets;r;p_{T};",cAdd.Data()),
 435                                      40,0.,2.,nBinPt,binLimitsPt);
 436         fHistList->Add(fh2PsiPt[ij][i]);
 437         fh2PhiPt[ij][i] = new TH2F(Form("fh2PhiPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;#phi;p_{T};",cAdd.Data()),
 438                                    nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 439         fHistList->Add(fh2PhiPt[ij][i]);
 440         fh2EtaPt[ij][i] = new TH2F(Form("fh2EtaPt%s_j%d",cAdd.Data(),i),Form("pt vs eta %s;#eta;p_{T};",cAdd.Data()),
 441                                    50,-1.,1.,nBinPt,binLimitsPt);
 442         fHistList->Add(fh2EtaPt[ij][i]);
 443         fh2AreaPt[ij][i] = new TH2F(Form("fh2AreaPt%s_j%d",cAdd.Data(),i),
 444                                     Form("pt vs area %s;area;p_{T};",
 445                                          cAdd.Data()),
 446                                     50,0.,1.,nBinPt,binLimitsPt);
 447         fHistList->Add(fh2AreaPt[ij][i]);
 448         fh2EtaArea[ij][i] = new TH2F(Form("fh2EtaArea%s_j%d",cAdd.Data(),i),
 449                                      Form("area vs eta %s;#eta;area;",
 450                                           cAdd.Data()),
 451                                      50,-1.,1.,50,0,1.);
 452         fHistList->Add(fh2EtaArea[ij][i]);
 453
 454         fh2PhiEta[ij][i] =  new TH2F(Form("fh2PhiEta%s_j%d",cAdd.Data(),i),Form("phi vs eta %s ;phi;p_{T};",cAdd.Data()),
 455                                      nBinPhi,binLimitsPhi,20,-1.,1.);
 456         fHistList->Add(fh2PhiEta[ij][i]);
 457
 458         fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
 459                                            Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
 460                                                 cAdd.Data()),
 461                                            200,0.,200.,nBinPt,binLimitsPt);
 462         fHistList->Add(fh2LTrackPtJetPt[ij][i]);
 463       }
 464
 465
 466       fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 467       fHistList->Add(fh1DijetMinv[ij]);
 468
 469       fh2DijetDeltaPhiPt[ij]       = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
 470       fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 471
 472       fh2DijetAsymPt[ij]            = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 473       fHistList->Add(fh2DijetAsymPt[ij]);
 474
 475       fh2DijetPt2vsPt1[ij]          = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
 476       fHistList->Add(fh2DijetPt2vsPt1[ij]);
 477       fh2DijetDifvsSum[ij]         = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
 478       fHistList->Add( fh2DijetDifvsSum[ij]);
 479     }
 480   // =========== Switch on Sumw2 for all histos ===========
 481   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 482     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 483     if (h1){
 484       h1->Sumw2();
 485       continue;
 486     }
 487     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 488     if(hn)hn->Sumw2();
 489   }
 490   TH1::AddDirectory(oldStatus);
 491 }
 492
 493 void AliAnalysisTaskJetSpectrum2::Init()
 494 {
 495   //
 496   // Initialization
 497   //
 498
 499   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 500
 501 }
 502
 503 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 504
 505
 506   Bool_t selected = kTRUE;
 507
 508   if(fUseGlobalSelection&&fEventSelectionMask==0){
 509     selected = AliAnalysisHelperJetTasks::Selected();
 510   }
 511   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 512     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 513   }
 514
 515   if(fEventClass>0){
 516     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 517   }
 518
 519   if(!selected){
 520     // no selection by the service task, we continue
 521     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 522     PostData(1, fHistList);
 523     return;
 524   }
 525
 526
 527   static AliAODEvent* aod = 0;
 528
 529   // take all other information from the aod we take the tracks from
 530   if(!aod){
 531    if(fUseAODTrackInput)aod = fAODIn;
 532    else aod = fAODOut;
 533   }
 534
 535
 536   //
 537   // Execute analysis for current event
 538   //
 539   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 540   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 541   if(!aodH){
 542     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 543     return;
 544   }
 545
 546   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 547   TClonesArray *aodRecJets = 0;
 548
 549   if(fAODOut&&!aodRecJets){
 550     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 551   }
 552   if(fAODExtension&&!aodRecJets){
 553     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 554   }
 555   if(fAODIn&&!aodRecJets){
 556     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 557   }
 558
 559
 560
 561   if(!aodRecJets){
 562     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 563     return;
 564   }
 565
 566   TClonesArray *aodGenJets = 0;
 567   if(fBranchGen.Length()>0){
 568     if(fAODOut&&!aodGenJets){
 569       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 570     }
 571     if(fAODExtension&&!aodGenJets){
 572       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 573     }
 574     if(fAODIn&&!aodGenJets){
 575       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 576     }
 577
 578     if(!aodGenJets){
 579       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 580       return;
 581     }
 582   }
 583
 584   TClonesArray *aodBackRecJets = 0;
 585   if(fBranchBkgRec.Length()>0){
 586     if(fAODOut&&!aodBackRecJets){
 587       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
 588     }
 589     if(fAODExtension&&!aodBackRecJets){
 590       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
 591     }
 592     if(fAODIn&&!aodBackRecJets){
 593       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
 594     }
 595
 596     if(!aodBackRecJets){
 597       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
 598       return;
 599     }
 600   }
 601
 602
 603   TClonesArray *aodBackGenJets = 0;
 604
 605   if(fBranchBkgGen.Length()>0){
 606     if(fAODOut&&!aodBackGenJets){
 607       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
 608     }
 609     if(fAODExtension&&!aodBackGenJets){
 610       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
 611     }
 612     if(fAODIn&&!aodBackGenJets){
 613       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
 614     }
 615
 616     if(!aodBackGenJets){
 617       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
 618       return;
 619     }
 620   }
 621
 622
 623   // new Scheme
 624   // first fill all the pure  histograms, i.e. full jets
 625   // in case of the correaltion limit to the n-leading jets
 626
 627   // reconstructed
 628
 629
 630   // generated
 631
 632
 633   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 634
 635
 636   TList genJetsList;         // full acceptance
 637   TList genJetsListCut;      // acceptance cut
 638   TList recJetsList;         // full acceptance
 639   TList recJetsListCut;      // acceptance cut
 640
 641
 642   GetListOfJets(&recJetsList,aodRecJets,0);
 643   GetListOfJets(&recJetsListCut,aodRecJets,1);
 644
 645   if(fBranchGen.Length()>0){
 646     GetListOfJets(&genJetsList,aodGenJets,0);
 647     GetListOfJets(&genJetsListCut,aodGenJets,1);
 648   }
 649
 650   Double_t eventW = 1;
 651   Double_t ptHard = 0;
 652   Double_t nTrials = 1; // Trials for MC trigger
 653   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 654
 655   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 656     // this is the part we only use when we have MC information
 657     AliMCEvent* mcEvent = MCEvent();
 658     //    AliStack *pStack = 0;
 659     if(!mcEvent){
 660       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 661       return;
 662     }
 663     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 664     if(pythiaGenHeader){
 665       nTrials = pythiaGenHeader->Trials();
 666       ptHard  = pythiaGenHeader->GetPtHard();
 667       int iProcessType = pythiaGenHeader->ProcessType();
 668       // 11 f+f -> f+f
 669       // 12 f+barf -> f+barf
 670       // 13 f+barf -> g+g
 671       // 28 f+g -> f+g
 672       // 53 g+g -> f+barf
 673       // 68 g+g -> g+g
 674       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 675       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 676     }
 677   }// (fAnalysisType&kMCESD)==kMCESD)
 678
 679
 680   // we simply fetch the tracks/mc particles as a list of AliVParticles
 681
 682   TList recParticles;
 683   TList genParticles;
 684
 685   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 686   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 687   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 688   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 689
 690   // Event control and counting ...
 691   // MC
 692   fh1PtHard->Fill(ptHard,eventW);
 693   fh1PtHardNoW->Fill(ptHard,1);
 694   fh1PtHardTrials->Fill(ptHard,nTrials);
 695
 696   // Real
 697   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 698     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 699   }
 700
 701
 702   Int_t recMult1 = recParticles.GetEntries();
 703   Int_t genMult1 = genParticles.GetEntries();
 704
 705   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
 706   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
 707
 708   fh2MultRec->Fill(recMult1,recMult2);
 709   fh2MultGen->Fill(genMult1,genMult2);
 710   fMultRec = recMult1;
 711   if(fMultRec<=0)fMultRec = recMult2;
 712   fMultGen = genMult1;
 713   if(fMultGen<=0)fMultGen = genMult2;
 714
 715   // the loops for rec and gen should be indentical... pass it to a separate
 716   // function ...
 717   // Jet Loop
 718   // Track Loop
 719   // Jet Jet Loop
 720   // Jet Track loop
 721
 722   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 723   FillJetHistos(recJetsList,recParticles,kJetRecFull);
 724   FillTrackHistos(recParticles,kJetRec);
 725
 726   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 727   FillJetHistos(genJetsList,genParticles,kJetGenFull);
 728   FillTrackHistos(genParticles,kJetGen);
 729
 730   // Here follows the jet matching part
 731   // delegate to a separated method?
 732
 733   FillMatchHistos(recJetsList,genJetsList);
 734
 735   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 736   PostData(1, fHistList);
 737 }
 738
 739 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 740
 741   if(iType>=kJetTypes){
 742     return;
 743   }
 744
 745   Int_t refMult = fMultRec;
 746   if(iType==kJetGen||iType==kJetGenFull){
 747     refMult = fMultGen;
 748   }
 749
 750   Int_t nJets = jetsList.GetEntries();
 751   fh1NJets[iType]->Fill(nJets);
 752
 753   if(nJets<=0)return;
 754
 755   Float_t ptOld = 1.E+32;
 756   Float_t pT0 = 0;
 757   Float_t pT1 = 0;
 758   Float_t phi0 = 0;
 759   Float_t phi1 = 0;
 760   Int_t ij0 = -1;
 761   Int_t ij1 = -1;
 762
 763
 764   for(int ij = 0;ij < nJets;ij++){
 765     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 766     Float_t ptJet = jet->Pt();
 767     fh1PtJetsIn[iType]->Fill(ptJet);
 768     fh2MultJetPt[iType]->Fill(refMult,ptJet);
 769     if(ptJet>ptOld){
 770       Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
 771     }
 772     ptOld = ptJet;
 773
 774     // find the dijets assume sorting and acceptance cut...
 775     if(ij==0){
 776       ij0 = ij;
 777       pT0 = ptJet;
 778       phi0 = jet->Phi();
 779       if(phi0<0)phi0+=TMath::Pi()*2.;
 780     }
 781     else if(ptJet>pT1){
 782       // take only the backward hemisphere??
 783       phi1 = jet->Phi();
 784       if(phi1<0)phi1+=TMath::Pi()*2.;
 785       Float_t dPhi = phi1 - phi0;
 786       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 787       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 788       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
 789         ij1 = ij;
 790         pT1 = ptJet;
 791       }
 792     }
 793     // fill jet histos for kmax jets
 794
 795       Float_t phiJet = jet->Phi();
 796       Float_t etaJet = jet->Eta();
 797       if(phiJet<0)phiJet+=TMath::Pi()*2.;
 798       fh1TmpRho->Reset();
 799       if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
 800       fh1PtIn[iType][kMaxJets]->Fill(ptJet);
 801       // fill leading jets...
 802       AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
 803       //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
 804       if(ptJet>10){
 805         if(ij<kMaxJets){
 806           fh2PhiEta[iType][ij]->Fill(phiJet,etaJet);
 807           fh2AreaPt[iType][ij]->Fill(jet->EffectiveAreaCharged(),ptJet);
 808           fh2EtaArea[iType][ij]->Fill(etaJet,jet->EffectiveAreaCharged());
 809         }
 810         fh2PhiEta[iType][kMaxJets]->Fill(phiJet,etaJet);
 811         fh2AreaPt[iType][kMaxJets]->Fill(jet->EffectiveAreaCharged(),ptJet);
 812         fh2EtaArea[iType][kMaxJets]->Fill(etaJet,jet->EffectiveAreaCharged());
 813       }
 814       if(ij<kMaxJets){
 815         fh2PhiPt[iType][ij]->Fill(phiJet,ptJet);
 816         fh2EtaPt[iType][ij]->Fill(etaJet,ptJet);
 817         if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
 818       }
 819       fh2PhiPt[iType][kMaxJets]->Fill(phiJet,ptJet);
 820       fh2EtaPt[iType][kMaxJets]->Fill(etaJet,ptJet);
 821       if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
 822
 823       if(particlesList.GetSize()&&ij<kMaxJets){
 824         // Particles... correlated with jets...
 825         for(int it = 0;it<particlesList.GetEntries();++it){
 826           AliVParticle *part = (AliVParticle*)particlesList.At(it);
 827           Float_t deltaR = jet->DeltaR(part);
 828           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
 829         }
 830         // fill the jet shapes
 831         Float_t rhoSum = 0;
 832         for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
 833           Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
 834           Float_t rho = fh1TmpRho->GetBinContent(ibx);
 835           rhoSum += rho;
 836           fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
 837           fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
 838         }
 839       }// if we have particles
 840   }// Jet Loop
 841
 842
 843   // do something with dijets...
 844   if(ij0>=0&&ij1>0){
 845     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
 846     Double_t ptJet0 = jet0->Pt();
 847     Double_t phiJet0 = jet0->Phi();
 848     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
 849
 850     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
 851     Double_t ptJet1 = jet1->Pt();
 852     Double_t phiJet1 = jet1->Phi();
 853     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
 854
 855     Float_t deltaPhi = phiJet0 - phiJet1;
 856     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
 857     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
 858     deltaPhi = TMath::Abs(deltaPhi);
 859     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
 860
 861     Float_t asym = 9999;
 862     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
 863       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
 864       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
 865       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
 866       Float_t minv = 2.*(jet0->P()*jet1->P()-
 867                          jet0->Px()*jet1->Px()-
 868                          jet0->Py()*jet1->Py()-
 869                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
 870       if(minv<0)minv=0; // prevent numerical instabilities
 871       minv = TMath::Sqrt(minv);
 872       fh1DijetMinv[iType]->Fill(minv);
 873   }
 874
 875
 876
 877   // count down the jets above thrueshold
 878   Int_t nOver = nJets;
 879   if(nOver>0){
 880     TIterator *jetIter = jetsList.MakeIterator();
 881     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
 882     if(tmpJet){
 883       Float_t pt = tmpJet->Pt();
 884       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
 885         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
 886         while(pt<ptCut&&tmpJet){
 887           nOver--;
 888           tmpJet = (AliAODJet*)(jetIter->Next());
 889           if(tmpJet){
 890             pt = tmpJet->Pt();
 891           }
 892         }
 893         if(nOver<=0)break;
 894         fh2NJetsPt[iType]->Fill(ptCut,nOver);
 895       }
 896     }
 897     delete jetIter;
 898   }
 899 }
 900
 901 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
 902
 903   Int_t nTrackOver = particlesList.GetSize();
 904   // do the same for tracks and jets
 905   if(nTrackOver>0){
 906     TIterator *trackIter = particlesList.MakeIterator();
 907     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
 908     Float_t pt = tmpTrack->Pt();
 909     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
 910       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
 911       while(pt<ptCut&&tmpTrack){
 912         nTrackOver--;
 913         tmpTrack = (AliVParticle*)(trackIter->Next());
 914         if(tmpTrack){
 915           pt = tmpTrack->Pt();
 916         }
 917       }
 918       if(nTrackOver<=0)break;
 919       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
 920     }
 921
 922
 923     trackIter->Reset();
 924     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
 925     Float_t sumPt = 0;
 926
 927     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
 928       Float_t tmpPt = tmpTrack->Pt();
 929       fh1PtTracksIn[iType]->Fill(tmpPt);
 930       sumPt += tmpPt;
 931       Float_t tmpPhi = tmpTrack->Phi();
 932       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 933       if(tmpTrack==leading){
 934         fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
 935         fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
 936         continue;
 937       }
 938     }
 939     fh1SumPtTrack[iType]->Fill(sumPt);
 940     delete trackIter;
 941   }
 942
 943 }
 944
 945
 946 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
 947
 948
 949   // Fill al the matching histos
 950   // It is important that the acceptances for the mathing are as large as possible
 951   // to avoid false matches on the edge of acceptance
 952   // therefore we add some extra matching jets as overhead
 953
 954   static TArrayI aGenIndex(recJetsList.GetEntries());
 955   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
 956
 957   static TArrayI aRecIndex(genJetsList.GetEntries());
 958   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
 959
 960   if(fDebug){
 961     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
 962     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
 963   }
 964   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
 965                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
 966                                             aGenIndex,aRecIndex,fDebug);
 967
 968   if(fDebug){
 969     for(int i = 0;i< aGenIndex.GetSize();++i){
 970       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
 971     }
 972     for(int i = 0;i< aRecIndex.GetSize();++i){
 973       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
 974     }
 975   }
 976
 977   Double_t container[6];
 978   Double_t containerPhiZ[6];
 979
 980   // loop over generated jets
 981   // consider the
 982   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
 983     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
 984     Double_t ptGen = genJet->Pt();
 985     Double_t phiGen = genJet->Phi();
 986     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 987     Double_t etaGen = genJet->Eta();
 988     container[3] = ptGen;
 989     container[4] = etaGen;
 990     container[5] = phiGen;
 991     fhnJetContainer[kStep0]->Fill(&container[3]);
 992     if(JetSelected(genJet)){
 993       fhnJetContainer[kStep1]->Fill(&container[3]);
 994       Int_t ir = aRecIndex[ig];
 995       if(ir>=0&&ir<recJetsList.GetEntries()){
 996         fhnJetContainer[kStep2]->Fill(&container[3]);
 997         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
 998         if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
 999         if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1000       }
1001     }
1002   }// loop over generated jets used for matching...
1003
1004
1005
1006   // loop over reconstructed jets
1007   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1008     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1009     Double_t etaRec = recJet->Eta();
1010     Double_t ptRec = recJet->Pt();
1011     Double_t phiRec = recJet->Phi();
1012     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1013     // do something with dijets...
1014
1015     container[0] = ptRec;
1016     container[1] = etaRec;
1017     container[2] = phiRec;
1018     containerPhiZ[0] = ptRec;
1019     containerPhiZ[1] = phiRec;
1020
1021     fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1022     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1023
1024     if(JetSelected(recJet)){
1025       fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1026       // Fill Correlation
1027       Int_t ig = aGenIndex[ir];
1028       if(ig>=0 && ig<genJetsList.GetEntries()){
1029         fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1030         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1031         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1032         Double_t ptGen  = genJet->Pt();
1033         Double_t phiGen = genJet->Phi();
1034         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1035         Double_t etaGen = genJet->Eta();
1036
1037         container[3] = ptGen;
1038         container[4] = etaGen;
1039         container[5] = phiGen;
1040         containerPhiZ[3] = ptGen;
1041         //
1042         // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1043         //
1044         if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1045         fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1046         fhnCorrelation->Fill(container);
1047         if(ptGen>0){
1048           Float_t delta = (ptRec-ptGen)/ptGen;
1049           fh2RelPtFGen->Fill(ptGen,delta);
1050           fh2PtFGen->Fill(ptGen,ptRec);
1051         }
1052         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1053       }
1054       else{
1055         containerPhiZ[3] = 0;
1056         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1057       }
1058     }// loop over reconstructed jets
1059   }
1060   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1061 }
1062
1063
1064 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1065   //
1066   // Create the particle container for the correction framework manager and
1067   // link it
1068   //
1069   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1070   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1071   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1072   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1073   const Double_t kZmin = 0., kZmax = 1;
1074
1075   // can we neglect migration in eta and phi?
1076   // phi should be no problem since we cover full phi and are phi symmetric
1077   // eta migration is more difficult  due to needed acceptance correction
1078   // in limited eta range
1079
1080   //arrays for the number of bins in each dimension
1081   Int_t iBin[kNvar];
1082   iBin[0] = 320; //bins in pt
1083   iBin[1] =  1; //bins in eta
1084   iBin[2] = 1; // bins in phi
1085
1086   //arrays for lower bounds :
1087   Double_t* binEdges[kNvar];
1088   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1089     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1090
1091   //values for bin lower bounds
1092   //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
1093   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1094   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1095   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1096
1097
1098   for(int i = 0;i < kMaxStep*2;++i){
1099     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1100     for (int k=0; k<kNvar; k++) {
1101       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1102     }
1103   }
1104   //create correlation matrix for unfolding
1105   Int_t thnDim[2*kNvar];
1106   for (int k=0; k<kNvar; k++) {
1107     //first half  : reconstructed
1108     //second half : MC
1109     thnDim[k]      = iBin[k];
1110     thnDim[k+kNvar] = iBin[k];
1111   }
1112
1113   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1114   for (int k=0; k<kNvar; k++) {
1115     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1116     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1117   }
1118   fhnCorrelation->Sumw2();
1119
1120   // for second correlation histogram
1121
1122
1123   const Int_t kNvarPhiZ   = 4;
1124   //arrays for the number of bins in each dimension
1125   Int_t iBinPhiZ[kNvarPhiZ];
1126   iBinPhiZ[0] = 80; //bins in pt
1127   iBinPhiZ[1] = 72; //bins in phi
1128   iBinPhiZ[2] = 20; // bins in Z
1129   iBinPhiZ[3] = 80; //bins in ptgen
1130
1131   //arrays for lower bounds :
1132   Double_t* binEdgesPhiZ[kNvarPhiZ];
1133   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1134     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1135
1136   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1137   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1138   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1139   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1140
1141   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1142   for (int k=0; k<kNvarPhiZ; k++) {
1143     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1144   }
1145   fhnCorrelationPhiZRec->Sumw2();
1146
1147
1148   // Add a histogram for Fake jets
1149
1150   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1151     delete [] binEdges[ivar];
1152
1153   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1154     delete [] binEdgesPhiZ[ivar];
1155
1156 }
1157
1158 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1159 {
1160 // Terminate analysis
1161 //
1162     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1163 }
1164
1165
1166 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1167
1168   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1169
1170   Int_t iCount = 0;
1171   if(type==kTrackAOD){
1172     AliAODEvent *aod = 0;
1173     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1174     else aod = AODEvent();
1175     if(!aod){
1176       return iCount;
1177     }
1178     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1179       AliAODTrack *tr = aod->GetTrack(it);
1180       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1181       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1182       if(tr->Pt()<fMinTrackPt)continue;
1183       if(fDebug>0){
1184         if(tr->Pt()>20){
1185           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1186           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1187           tr->Print();
1188           //    tr->Dump();
1189           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1190           if(fESD){
1191             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1192             esdTr->Print("");
1193             // esdTr->Dump();
1194           }
1195         }
1196       }
1197       list->Add(tr);
1198       iCount++;
1199     }
1200   }
1201   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1202     AliMCEvent* mcEvent = MCEvent();
1203     if(!mcEvent)return iCount;
1204     // we want to have alivpartilces so use get track
1205     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1206       if(!mcEvent->IsPhysicalPrimary(it))continue;
1207       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1208       if(part->Pt()<fMinTrackPt)continue;
1209       if(type == kTrackKineAll){
1210         list->Add(part);
1211         iCount++;
1212       }
1213       else if(type == kTrackKineCharged){
1214         if(part->Particle()->GetPDG()->Charge()==0)continue;
1215         list->Add(part);
1216         iCount++;
1217       }
1218     }
1219   }
1220   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1221     AliAODEvent *aod = 0;
1222     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1223     else aod = AODEvent();
1224     if(!aod)return iCount;
1225     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1226     if(!tca)return iCount;
1227     for(int it = 0;it < tca->GetEntriesFast();++it){
1228       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1229       if(!part)continue;
1230       if(part->Pt()<fMinTrackPt)continue;
1231       if(!part->IsPhysicalPrimary())continue;
1232       if(type == kTrackAODMCAll){
1233         list->Add(part);
1234         iCount++;
1235       }
1236       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1237         if(part->Charge()==0)continue;
1238         if(kTrackAODMCCharged){
1239           list->Add(part);
1240         }
1241         else {
1242           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1243           list->Add(part);
1244         }
1245         iCount++;
1246       }
1247     }
1248   }// AODMCparticle
1249   list->Sort();
1250   return iCount;
1251
1252 }
1253
1254
1255
1256 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1257   Bool_t selected = false;
1258
1259   if(!jet)return selected;
1260
1261   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1262     selected = kTRUE;
1263   }
1264   return selected;
1265
1266 }
1267
1268 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1269
1270   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1271   Int_t iCount = 0;
1272
1273   if(!jarray){
1274     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1275     return iCount;
1276   }
1277
1278
1279   for(int ij=0;ij<jarray->GetEntries();ij++){
1280     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1281     if(!jet)continue;
1282     if(type==0){
1283       // No cut at all, main purpose here is sorting
1284       list->Add(jet);
1285       iCount++;
1286     }
1287     else if(type == 1){
1288       // eta cut
1289       if(JetSelected(jet)){
1290         list->Add(jet);
1291         iCount++;
1292       }
1293     }
1294   }
1295
1296   list->Sort();
1297   return iCount;
1298
1299 }
1300
1301
1302 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1303   if(!jets)return 0;
1304
1305   Int_t refMult = 0;
1306   for(int ij = 0;ij < jets->GetEntries();++ij){
1307     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1308     if(!jet)continue;
1309     TRefArray *refs = jet->GetRefTracks();
1310     if(!refs)continue;
1311     refMult += refs->GetEntries();
1312   }
1313   return refMult;
1314
1315 }
1316
1317
1318 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1319   if(!jet)return 0;
1320   TRefArray *refs = jet->GetRefTracks();
1321   if(!refs) return 0;
1322   AliVParticle *leading = 0;
1323   Float_t fMaxPt = 0;
1324   for(int i = 0;i<refs->GetEntries();i++){
1325     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1326     if(!tmp)continue;
1327     if(tmp->Pt()>fMaxPt){
1328       leading = tmp;
1329       fMaxPt = tmp->Pt();
1330     }
1331   }
1332   return leading;
1333 }
1334
1335
1336 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1337   if(!jet)return 0;
1338   if(!list) return 0;
1339   AliVParticle *leading = 0;
1340   Float_t fMaxPt = 0;
1341   for(int i = 0;i<list->GetEntries();i++){
1342     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1343     if(!tmp)continue;
1344     if(jet->DeltaR(tmp)>r)continue;
1345     if(tmp->Pt()>fMaxPt){
1346       leading = tmp;
1347       fMaxPt = tmp->Pt();
1348     }
1349   }
1350   return leading;
1351 }