PWGJE/AliAnalysisTaskJetSpectrum2.cxx

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