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