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