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 ptOld = 1.E+32;
 988   Float_t pT0 = 0;
 989   Float_t pT1 = 0;
 990   Float_t phi0 = 0;
 991   Float_t phi1 = 0;
 992   Int_t ij0 = -1;
 993   Int_t ij1 = -1;
 994
 995   Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
 996   Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
 997
 998   var1[2] = fCentrality;
 999   var1b[1] = fCentrality;
1000   var1[3] = refMult;
1001
1002
1003
1004   Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
1005   var2[2] = fCentrality;
1006
1007   for(int ij = 0;ij < nJets;ij++){
1008     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
1009     Float_t ptJet = jet->Pt();
1010
1011
1012     if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
1013
1014     var1b[0] = ptJet;
1015     if(jet->Trigger()&fJetTriggerBestMask){
1016       fh1PtJetsInBest[iType]->Fill(ptJet);
1017       for(int it = 0;it <fNTrigger;it++){
1018         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1019           var1b[2] = it;
1020           fhnJetPtBest[iType]->Fill(var1b);
1021         }
1022       }
1023     }
1024     if(jet->Trigger()&fJetTriggerExcludeMask){
1025       fh1PtJetsInRej[iType]->Fill(ptJet);
1026       for(int it = 0;it <fNTrigger;it++){
1027         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1028           var1b[2] = it;
1029           fhnJetPtRej[iType]->Fill(var1b);
1030         }
1031       }
1032       continue;
1033     }
1034     fh1PtJetsIn[iType]->Fill(ptJet);
1035     ptOld = ptJet;
1036
1037     // find the dijets assume sorting and acceptance cut...
1038     if(ij==0){
1039       ij0 = ij;
1040       pT0 = ptJet;
1041       phi0 = jet->Phi();
1042       if(phi0<0)phi0+=TMath::Pi()*2.;
1043     }
1044     else if(ptJet>pT1){
1045       // take only the backward hemisphere??
1046       phi1 = jet->Phi();
1047       if(phi1<0)phi1+=TMath::Pi()*2.;
1048       Float_t dPhi = phi1 - phi0;
1049       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1050       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1051       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1052         ij1 = ij;
1053         pT1 = ptJet;
1054       }
1055     }
1056     // fill jet histos for kmax jets
1057     Float_t phiJet = jet->Phi();
1058     Float_t etaJet = jet->Eta();
1059     if(phiJet<0)phiJet+=TMath::Pi()*2.;
1060     fh1TmpRho->Reset();
1061     if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1062
1063     fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1064     // fill leading jets...
1065     AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1066     //AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1067     Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1068     Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1069     AliVParticle *tt = (AliVParticle*)particlesList.At(0);
1070     Float_t ttphi=tt->Phi();
1071     if(ttphi<0)ttphi+=TMath::Pi()*2.;
1072     Float_t ttpt=tt->Pt();
1073     Int_t phiBintt = GetPhiBin(ttphi-fRPAngle);
1074     Double_t dphitrigjet=RelativePhi(ttphi,phiJet);
1075     if(fTRP==1){
1076       if(TMath::Abs(dphitrigjet)<TMath::Pi()-0.6) continue;
1077       var1[1] = ptJet;
1078       var1[4] = phiBintt;
1079       var1[5] = jet->EffectiveAreaCharged();
1080       var1[7] = ttpt;
1081       var1[8] = CheckAcceptance(phiJet,etaJet);
1082     }
1083
1084     if(fTRP==0){
1085       var1[1] = ptJet;
1086       var1[4] = phiBin;
1087       var1[5] = jet->EffectiveAreaCharged();
1088       var1[7] = ptLead;
1089       var1[8] = CheckAcceptance(phiJet,etaJet);
1090     }
1091
1092     //jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading
1093     var2[1] = ptJet;
1094     var2[3] = etaJet;
1095     var2[4] = phiJet;
1096     var2[5] = jet->EffectiveAreaCharged();
1097     var2[7] = var1[8];
1098     var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1099
1100     if(ij<kMaxJets){
1101       fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1102       var1[0] = ij;
1103       var2[0] = ij;
1104       for(int it = 0;it <fNTrigger;it++){
1105         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1106           var1[6] = it;
1107           var2[6] = it;
1108           fhnJetPt[iType]->Fill(var1);
1109           fhnJetPtQA[iType]->Fill(var2);
1110         }
1111       }
1112     }
1113     var1[0] = kMaxJets;// fill for all jets
1114     var2[0] = kMaxJets;// fill for all jets
1115     for(int it = 0;it <fNTrigger;it++){
1116       if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1117         var1[6] = it;
1118         var2[6] = it;
1119         fhnJetPt[iType]->Fill(var1);
1120         fhnJetPtQA[iType]->Fill(var2);
1121       }
1122     }
1123
1124     fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1125
1126     if(particlesList.GetSize()&&ij<kMaxJets){
1127       // Particles... correlated with jets...
1128       for(int it = 0;it<particlesList.GetEntries();++it){
1129         AliVParticle *part = (AliVParticle*)particlesList.At(it);
1130         Float_t deltaR = jet->DeltaR(part);
1131         if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1132       }
1133       // fill the jet shapes
1134     }// if we have particles
1135   }// Jet Loop
1136
1137
1138   // do something with dijets...
1139   if(ij0>=0&&ij1>0){
1140     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1141     Double_t ptJet0 = jet0->Pt();
1142     Double_t phiJet0 = jet0->Phi();
1143     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1144
1145     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1146     Double_t ptJet1 = jet1->Pt();
1147     Double_t phiJet1 = jet1->Phi();
1148     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1149
1150     Float_t deltaPhi = phiJet0 - phiJet1;
1151     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1152     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1153     deltaPhi = TMath::Abs(deltaPhi);
1154     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1155
1156     Float_t asym = 9999;
1157     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1158     fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1159     fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1160     fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1161     Float_t minv = 2.*(jet0->P()*jet1->P()-
1162                        jet0->Px()*jet1->Px()-
1163                        jet0->Py()*jet1->Py()-
1164                        jet0->Pz()*jet1->Pz());    // assume mass == 0;
1165     if(minv<0)minv=0; // prevent numerical instabilities
1166     minv = TMath::Sqrt(minv);
1167     fh1DijetMinv[iType]->Fill(minv);
1168   }
1169
1170
1171
1172   // count down the jets above thrueshold
1173   Int_t nOver = nJets;
1174   if(nOver>0){
1175     TIterator *jetIter = jetsList.MakeIterator();
1176     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1177     if(tmpJet){
1178       Float_t pt = tmpJet->Pt();
1179       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1180         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1181         while(pt<ptCut&&tmpJet){
1182           nOver--;
1183           tmpJet = (AliAODJet*)(jetIter->Next());
1184           if(tmpJet){
1185             pt = tmpJet->Pt();
1186           }
1187         }
1188         if(nOver<=0)break;
1189         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1190       }
1191     }
1192     delete jetIter;
1193   }
1194 }
1195
1196 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1197
1198   if(fFlagJetType[iType]<=0)return;
1199   Int_t refMult = fMultRec;
1200   if(iType==kJetGen||iType==kJetGenFull){
1201     refMult = fMultGen;
1202
1203   }
1204
1205   //
1206   Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge
1207   var3[2] = fCentrality;
1208   var3[3] = refMult;
1209   Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1210   var4[2] = fCentrality;
1211   Int_t nTrackOver = particlesList.GetSize();
1212   // do the same for tracks and jets
1213   if(nTrackOver>0){
1214     TIterator *trackIter = particlesList.MakeIterator();
1215     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1216     Float_t pt = tmpTrack->Pt();
1217     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1218       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1219       while(pt<ptCut&&tmpTrack){
1220         nTrackOver--;
1221         tmpTrack = (AliVParticle*)(trackIter->Next());
1222         if(tmpTrack){
1223           pt = tmpTrack->Pt();
1224         }
1225       }
1226       if(nTrackOver<=0)break;
1227       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1228     }
1229
1230
1231     trackIter->Reset();
1232     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1233     Float_t sumPt = 0;
1234
1235     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1236       Float_t tmpPt = tmpTrack->Pt();
1237       fh1PtTracksIn[iType]->Fill(tmpPt);
1238       fh1PtTracksInLow[iType]->Fill(tmpPt);
1239
1240       sumPt += tmpPt;
1241       Float_t tmpPhi = tmpTrack->Phi();
1242       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1243
1244
1245       Float_t phiRP = tmpPhi-fRPAngle;
1246       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1247       if(phiRP<0)phiRP += TMath::Pi();
1248       if(phiRP<0)phiRP += TMath::Pi();
1249       const float allPhi = -1./180.*TMath::Pi();
1250
1251       if(tmpPt<100){
1252         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1253         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1254
1255         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1256         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1257       }
1258       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1259       var3[0] = 1;
1260       var3[1] = tmpPt;
1261       var3[4] = phiBin;
1262       var3[6] = tmpTrack->Charge();
1263
1264       var4[0] = 1;
1265       var4[1] = tmpPt;
1266       var4[3] = tmpTrack->Eta();
1267       var4[4] = tmpPhi;
1268       var4[5] = tmpTrack->Charge();
1269
1270
1271       for(int it = 0;it <fNTrigger;it++){
1272         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1273           var3[0] = 1;
1274           var3[5] = it;
1275           fhnTrackPt[iType]->Fill(var3);
1276           if(tmpTrack==leading){
1277             var3[0] = 0;
1278             fhnTrackPt[iType]->Fill(var3);
1279           }
1280         }
1281       }
1282
1283
1284
1285
1286       fhnTrackPtQA[iType]->Fill(var4);
1287       if(tmpTrack==leading){
1288         var4[0] = 0;
1289         fhnTrackPtQA[iType]->Fill(var4);
1290         continue;
1291       }
1292
1293
1294
1295
1296
1297     }
1298     fh1SumPtTrack[iType]->Fill(sumPt);
1299     delete trackIter;
1300   }
1301
1302 }
1303
1304
1305 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1306
1307
1308   // Fill al the matching histos
1309   // It is important that the acceptances for the mathing are as large as possible
1310   // to avoid false matches on the edge of acceptance
1311   // therefore we add some extra matching jets as overhead
1312
1313   static TArrayI aGenIndex(recJetsList.GetEntries());
1314   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1315
1316   static TArrayI aRecIndex(genJetsList.GetEntries());
1317   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1318
1319   if(fDebug){
1320     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1321     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1322   }
1323   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1324                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1325                                             aGenIndex,aRecIndex,fDebug);
1326
1327   if(fDebug){
1328     for(int i = 0;i< aGenIndex.GetSize();++i){
1329       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1330     }
1331     for(int i = 0;i< aRecIndex.GetSize();++i){
1332       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1333     }
1334   }
1335
1336   Double_t container[8];
1337   Double_t containerRec[4];
1338   Double_t containerGen[4];
1339
1340   // loop over generated jets
1341   // consider the
1342   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1343     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1344     Double_t ptGen = genJet->Pt();
1345     Double_t phiGen = genJet->Phi();
1346     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1347     Double_t etaGen = genJet->Eta();
1348     containerGen[0] = ptGen;
1349     containerGen[1] = etaGen;
1350     containerGen[2] = phiGen;
1351     containerGen[3] = genJet->GetPtLeading();
1352
1353     fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1354
1355     if(JetSelected(genJet))
1356       fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1357
1358     Int_t ir = aRecIndex[ig];
1359     if(ir>=0&&ir<recJetsList.GetEntries()){
1360       AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1361
1362       fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1363
1364       if(JetSelected(genJet)){
1365         fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1366
1367         if(JetSelected(recJet)) {
1368           fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1369
1370         }
1371         containerGen[3] = recJet->GetPtLeading();
1372         fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level
1373       }
1374     }
1375   }// loop over generated jets used for matching...
1376
1377
1378
1379   // loop over reconstructed jets
1380   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1381     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1382     Double_t etaRec = recJet->Eta();
1383     Double_t ptRec = recJet->Pt();
1384     Double_t phiRec = recJet->Phi();
1385     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1386
1387     containerRec[0] = ptRec;
1388     containerRec[1] = etaRec;
1389     containerRec[2] = phiRec;
1390     containerRec[3] = recJet->GetPtLeading();
1391
1392     container[0] = ptRec;
1393     container[1] = etaRec;
1394     container[2] = phiRec;
1395     container[3] = recJet->GetPtLeading();
1396
1397     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1398
1399     if(JetSelected(recJet)){
1400       fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1401       // Fill Correlation
1402       Int_t ig = aGenIndex[ir];
1403       if(ig>=0 && ig<genJetsList.GetEntries()){
1404         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1405         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1406         Double_t ptGen  = genJet->Pt();
1407         Double_t phiGen = genJet->Phi();
1408         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1409         Double_t etaGen = genJet->Eta();
1410
1411         containerGen[0] = ptGen;
1412         containerGen[1] = etaGen;
1413         containerGen[2] = phiGen;
1414         containerGen[3] = genJet->GetPtLeading();
1415
1416         container[4] = ptGen;
1417         container[5] = etaGen;
1418         container[6] = phiGen;
1419         container[7] = genJet->GetPtLeading();
1420
1421         fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1422
1423         fhnCorrelation->Fill(container);
1424
1425         Float_t delta = (ptRec-ptGen)/ptGen;
1426         fh2RelPtFGen->Fill(ptGen,delta);
1427         fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading());
1428         fh2PtFGen->Fill(ptGen,ptRec);
1429
1430         fh2deltaPt1Pt2->Fill(ptRec-ptGen,ptRec,ptGen);
1431
1432       }
1433     }// loop over reconstructed jets
1434   }
1435   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1436 }
1437
1438
1439 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1440   //
1441   // Create the particle container for the correction framework manager and
1442   // link it
1443   //
1444   const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
1445   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1446   const Double_t kEtamin = -1.0, kEtamax = 1.0;
1447   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1448   const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1449
1450   // can we neglect migration in eta and phi?
1451   // phi should be no problem since we cover full phi and are phi symmetric
1452   // eta migration is more difficult  due to needed acceptance correction
1453   // in limited eta range
1454
1455   //arrays for the number of bins in each dimension
1456   Int_t iBin[kNvar];
1457   iBin[0] = 125; //bins in pt
1458   iBin[1] =  4; //bins in eta
1459   iBin[2] =  1; // bins in phi
1460   iBin[3] =  10; // bins in leading track Pt
1461
1462   //arrays for lower bounds :
1463   Double_t* binEdges[kNvar];
1464   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1465     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1466
1467   //values for bin lower bounds
1468   //  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);
1469   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1470   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1471   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1472   binEdges[3][0]= kPtLeadingTrackPtMin;
1473   binEdges[3][1]= 1.;
1474   binEdges[3][2]= 2.;
1475   binEdges[3][3]= 3.;
1476   binEdges[3][4]= 4.;
1477   binEdges[3][5]= 5.;
1478   binEdges[3][6]= 6.;
1479   binEdges[3][7]= 8.;
1480   binEdges[3][8]= 10.;
1481   binEdges[3][9]= 12.;
1482   binEdges[3][10]= kPtLeadingTrackPtMax;
1483
1484   fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1485   for (int k=0; k<kNvar; k++) {
1486     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1487   }
1488
1489   //create correlation matrix for unfolding
1490   Int_t thnDim[2*kNvar];
1491   for (int k=0; k<kNvar; k++) {
1492     //first half  : reconstructed
1493     //second half : MC
1494     thnDim[k]      = iBin[k];
1495     thnDim[k+kNvar] = iBin[k];
1496   }
1497
1498   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1499   for (int k=0; k<kNvar; k++) {
1500     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1501     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1502   }
1503
1504   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1505     delete [] binEdges[ivar];
1506
1507
1508 }
1509
1510 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1511 {
1512   // Terminate analysis
1513   //
1514   if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1515 }
1516
1517
1518 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1519
1520   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1521
1522   Int_t iCount = 0;
1523   if(type==kTrackAOD){
1524     AliAODEvent *aod = 0;
1525     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1526     else aod = AODEvent();
1527     if(!aod){
1528       return iCount;
1529     }
1530     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1531       AliAODTrack *tr = aod->GetTrack(it);
1532       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1533       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1534       if(tr->Pt()<fMinTrackPt)continue;
1535       if(fDebug>0){
1536         if(tr->Pt()>20){
1537           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1538           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1539           tr->Print();
1540           //    tr->Dump();
1541           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1542           if(fESD){
1543             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1544             esdTr->Print("");
1545             // esdTr->Dump();
1546           }
1547         }
1548       }
1549       list->Add(tr);
1550       iCount++;
1551     }
1552   }
1553   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1554     AliMCEvent* mcEvent = MCEvent();
1555     if(!mcEvent)return iCount;
1556     // we want to have alivpartilces so use get track
1557     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1558       if(!mcEvent->IsPhysicalPrimary(it))continue;
1559       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1560       if(part->Pt()<fMinTrackPt)continue;
1561       if(type == kTrackKineAll){
1562         list->Add(part);
1563         iCount++;
1564       }
1565       else if(type == kTrackKineCharged){
1566         if(part->Particle()->GetPDG()->Charge()==0)continue;
1567         list->Add(part);
1568         iCount++;
1569       }
1570     }
1571   }
1572   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1573     AliAODEvent *aod = 0;
1574     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1575     else aod = AODEvent();
1576     if(!aod)return iCount;
1577     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1578     if(!tca)return iCount;
1579     for(int it = 0;it < tca->GetEntriesFast();++it){
1580       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1581       if(!part)continue;
1582       if(part->Pt()<fMinTrackPt)continue;
1583       if(!part->IsPhysicalPrimary())continue;
1584       if(type == kTrackAODMCAll){
1585         list->Add(part);
1586         iCount++;
1587       }
1588       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1589         if(part->Charge()==0)continue;
1590         if(kTrackAODMCCharged){
1591           list->Add(part);
1592         }
1593         else {
1594           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1595           list->Add(part);
1596         }
1597         iCount++;
1598       }
1599     }
1600   }// AODMCparticle
1601   list->Sort();
1602   return iCount;
1603
1604 }
1605
1606
1607 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1608   AliAODEvent *aod = 0;
1609   if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1610   else aod = AODEvent();
1611   if(!aod){
1612     return 101;
1613   }
1614   return aod->GetHeader()->GetCentrality();
1615 }
1616
1617
1618
1619 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1620   Bool_t selected = false;
1621
1622   if(!jet)return selected;
1623
1624   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1625     selected = kTRUE;
1626   }
1627   return selected;
1628
1629 }
1630
1631 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1632
1633   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1634   Int_t iCount = 0;
1635
1636   if(!jarray){
1637     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1638     return iCount;
1639   }
1640
1641
1642   for(int ij=0;ij<jarray->GetEntries();ij++){
1643     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1644     if(!jet)continue;
1645     if(type==0){
1646       // No cut at all, main purpose here is sorting
1647       list->Add(jet);
1648       iCount++;
1649     }
1650     else if(type == 1){
1651       // eta cut
1652       if(JetSelected(jet)){
1653         list->Add(jet);
1654         iCount++;
1655       }
1656     }
1657   }
1658
1659   list->Sort();
1660   return iCount;
1661
1662 }
1663
1664
1665 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1666   if(!jets)return 0;
1667
1668   Int_t refMult = 0;
1669   for(int ij = 0;ij < jets->GetEntries();++ij){
1670     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1671     if(!jet)continue;
1672     TRefArray *refs = jet->GetRefTracks();
1673     if(!refs)continue;
1674     refMult += refs->GetEntries();
1675   }
1676   return refMult;
1677
1678 }
1679
1680
1681 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1682   if(!jet)return 0;
1683   TRefArray *refs = jet->GetRefTracks();
1684   if(!refs) return 0;
1685   AliVParticle *leading = 0;
1686   Float_t fMaxPt = 0;
1687   for(int i = 0;i<refs->GetEntries();i++){
1688     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1689     if(!tmp)continue;
1690     if(tmp->Pt()>fMaxPt){
1691       leading = tmp;
1692       fMaxPt = tmp->Pt();
1693     }
1694   }
1695   return leading;
1696 }
1697
1698
1699 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1700   if(!jet)return 0;
1701   if(!list) return 0;
1702   AliVParticle *leading = 0;
1703   Float_t fMaxPt = 0;
1704   for(int i = 0;i<list->GetEntries();i++){
1705     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1706     if(!tmp)continue;
1707     if(jet->DeltaR(tmp)>r)continue;
1708     if(tmp->Pt()>fMaxPt){
1709       leading = tmp;
1710       fMaxPt = tmp->Pt();
1711     }
1712   }
1713   return leading;
1714 }
1715
1716
1717 Double_t AliAnalysisTaskJetSpectrum2::RelativePhi(Double_t mphi,Double_t vphi){
1718
1719   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1720   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1721   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1722   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1723   double dphi = mphi-vphi;
1724   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1725   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1726   return dphi;//dphi in [-Pi, Pi]
1727 }
1728
1729
1730
1731
1732
1733
1734 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1735 {
1736   Int_t phibin=-1;
1737   if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1738   Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1739   phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1740   if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1741   return phibin;
1742 }
1743
1744 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1745   //
1746   // set number of trigger bins
1747   //
1748   if(n>0){
1749     fNTrigger = n;
1750     delete [] fTriggerName;
1751     fTriggerName = new TString [fNTrigger];
1752     delete [] fTriggerBit;fTriggerBit = 0;
1753     fTriggerBit = new UInt_t [fNTrigger];
1754   }
1755   else{
1756     fNTrigger = 0;
1757   }
1758 }
1759
1760
1761 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1762   //
1763   // set trigger bin
1764   //
1765   if(i<fNTrigger){
1766     fTriggerBit[i] = it;
1767     fTriggerName[i] =  c;
1768   }
1769 }
1770
1771
1772 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1773   //
1774   // Set number of acceptance bins
1775   //
1776
1777   if(n>0){
1778     fNAcceptance = n;
1779     delete [] fAcceptancePhiMin;
1780     delete [] fAcceptancePhiMax;
1781     delete [] fAcceptanceEtaMin;
1782     delete [] fAcceptanceEtaMax;
1783
1784     fAcceptancePhiMin = new Float_t[fNAcceptance];
1785     fAcceptancePhiMax = new Float_t[fNAcceptance];
1786     fAcceptanceEtaMin = new Float_t[fNAcceptance];
1787     fAcceptanceEtaMax = new Float_t[fNAcceptance];
1788   }
1789   else{
1790     fNTrigger = 0;
1791   }
1792 }
1793
1794 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1795   //
1796   // Set acceptance bins
1797   //
1798
1799   if(i<fNAcceptance){
1800     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1801
1802     fAcceptancePhiMin[i] = phiMin;
1803     fAcceptancePhiMax[i] = phiMax;
1804     fAcceptanceEtaMin[i] = etaMin;
1805     fAcceptanceEtaMax[i] = etaMax;
1806   }
1807   else{
1808     fNTrigger = 0;
1809   }
1810 }
1811
1812
1813 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1814
1815   //
1816   // return aceptnace bin do no use ovelapping bins
1817   //
1818
1819   for(int i = 0;i<fNAcceptance;i++){
1820     if(phi<fAcceptancePhiMin[i])continue;
1821     if(phi>fAcceptancePhiMax[i])continue;
1822     if(eta<fAcceptanceEtaMin[i])continue;
1823     if(eta>fAcceptanceEtaMax[i])continue;
1824     return i;
1825   }
1826   // catch the rest
1827   return fNAcceptance;
1828 }
1829
1830 Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
1831
1832   // invert the correction
1833   AliAODJet *jet = (AliAODJet*)list.At(0);  // highest pt jet
1834   if(!jet)return -1;
1835   if(jet->EffectiveAreaCharged()<=0)return -1;
1836   Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1837   return rho;
1838
1839 }
1840
1841
1842
1843 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1844   //
1845   delete [] fTriggerBit;
1846   delete [] fTriggerName;
1847
1848   delete [] fAcceptancePhiMin;
1849   delete [] fAcceptancePhiMax;
1850   delete [] fAcceptanceEtaMin;
1851   delete [] fAcceptanceEtaMax;
1852
1853
1854
1855 }