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