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