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