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