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