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