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