PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include <Riostream.h>
  17
  18 #include <TFile.h>
  19 #include <TClonesArray.h>
  20 #include <TLorentzVector.h>
  21 #include <TProfile.h>
  22 #include <TList.h>
  23 #include <TH1.h>
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TH3F.h>
  27
  28 #include "AliAnalysisTaskFragmentationFunction.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliInputEventHandler.h"
  31 #include "AliAODHandler.h"
  32 #include "AliAODTrack.h"
  33 #include "AliJetHeader.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODJet.h"
  36 #include "AliAODDiJet.h"
  37 #include "AliGenPythiaEventHeader.h"
  38 #include "AliAnalysisHelperJetTasks.h"
  39 #include "AliMCEvent.h"
  40 #include "AliMCParticle.h"
  41 #include "AliAODMCParticle.h"
  42
  43 ClassImp(AliAnalysisTaskFragmentationFunction)
  44
  45 //#######################################################################
  46 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction():
  47   AliAnalysisTaskSE(),
  48   fJetHeaderRec(0x0),
  49   fJetHeaderGen(0x0),
  50   fAOD(0x0),
  51   fBranchRec(""),
  52   fBranchGen(""),
  53   fUseAODInput(kFALSE),
  54   fUseAODJetInput(kFALSE),
  55   fUseAODTrackInput(kFALSE),
  56   fUseAODMCInput(kFALSE),
  57   fUseGlobalSelection(kFALSE),
  58   fUseExternalWeightOnly(kFALSE),
  59   fLimitGenJetEta(kFALSE),
  60   fFilterMask(0),
  61   fAnalysisType(0),
  62   fTrackTypeRec(kTrackUndef),
  63   fTrackTypeGen(kTrackUndef),
  64   fAvgTrials(1.),
  65   fExternalWeight(1.),
  66   fRecEtaWindow(0.5),
  67   fR(0.),
  68   fdRdNdxi(0.7),
  69   fPartPtCut(0.),
  70   fEfactor(1.),
  71   fNff(5),
  72   fNim(5),
  73   fList(0x0),
  74   fGlobVar(1),
  75   // Number of energy bins
  76   fnEBin(6),
  77   fEmin(10.),
  78   fEmax(70.),
  79   fnEInterval(6),
  80   // Number of radius bins
  81   fnRBin(10),
  82   fRmin(0.1),
  83   fRmax(1.),
  84   fnRInterval(9),
  85   // eta histograms
  86   fnEtaHBin(50),
  87   fEtaHBinMin(-0.9),
  88   fEtaHBinMax(0.9),
  89   // phi histograms
  90   fnPhiHBin(60),
  91   fPhiHBinMin(0.),
  92   fPhiHBinMax(2*TMath::Pi()),
  93   // pt histograms
  94   fnPtHBin(300),
  95   fPtHBinMin(0.),
  96   fPtHBinMax(300.),
  97   // E histograms
  98   fnEHBin(300),
  99   fEHBinMin(0.),
 100   fEHBinMax(300.),
 101   // Xi histograms
 102   fnXiHBin(27),
 103   fXiHBinMin(0.),
 104   fXiHBinMax(9.),
 105   // Pthad histograms
 106   fnPthadHBin(60),
 107   fPthadHBinMin(0.),
 108   fPthadHBinMax(30.),
 109   // z histograms
 110   fnZHBin(30),
 111   fZHBinMin(0.),
 112   fZHBinMax(1.),
 113   // theta histograms
 114   fnThetaHBin(200),
 115   fThetaHBinMin(-0.5),
 116   fThetaHBinMax(1.5),
 117   fnCosThetaHBin(100),
 118   fcosThetaHBinMin(0.),
 119   fcosThetaHBinMax(1.),
 120   // kT histograms
 121   fnkTHBin(25),
 122   fkTHBinMin(0.),
 123   fkTHBinMax(5.),
 124   // kT histograms
 125   fnRHBin(10),
 126   fRHBinMin(0),
 127   fRHBinMax(1),
 128   // pt trig
 129   fnPtTrigBin(10),
 130   //Histograms
 131   fEtaMonoJet1H(0x0),
 132   fPhiMonoJet1H(0x0),
 133   fPtMonoJet1H(0x0),
 134   fEMonoJet1H(0x0),
 135   fdNdXiMonoJet1H(0x0),
 136   fdNdPtMonoJet1H(0x0),
 137   fdNdZMonoJet1H(0x0),
 138   fdNdThetaMonoJet1H(0x0),
 139   fdNdcosThetaMonoJet1H(0x0),
 140   fdNdkTMonoJet1H(0x0),
 141   fdNdpTvsZMonoJet1H(0x0),
 142   fShapeMonoJet1H(0x0),
 143   fNMonoJet1sH(0x0),
 144   fThetaPtPartMonoJet1H(0x0),
 145   fcosThetaPtPartMonoJet1H(0x0),
 146   fkTPtPartMonoJet1H(0x0),
 147   fThetaPtJetMonoJet1H(0x0),
 148   fcosThetaPtJetMonoJet1H(0x0),
 149   fkTPtJetMonoJet1H(0x0),
 150   fpTPtJetMonoJet1H(0x0),
 151   farrayEmin(0x0),
 152   farrayEmax(0x0),
 153   farrayRadii(0x0),
 154   farrayPtTrigmin(0x0),
 155   farrayPtTrigmax(0x0),
 156   // Track control plots
 157   fptAllTracks(0x0),
 158   fetaAllTracks(0x0),
 159   fphiAllTracks(0x0),
 160   fetaphiptAllTracks(0x0),
 161   fetaphiAllTracks(0x0),
 162   fptAllTracksCut(0x0),
 163   fetaAllTracksCut(0x0),
 164   fphiAllTracksCut(0x0),
 165   fetaphiptAllTracksCut(0x0),
 166   fetaphiAllTracksCut(0x0),
 167   fptTracks(0x0),
 168   fetaTracks(0x0),
 169   fphiTracks(0x0),
 170   fdetaTracks(0x0),
 171   fdphiTracks(0x0),
 172   fetaphiptTracks(0x0),
 173   fetaphiTracks(0x0),
 174   fdetadphiTracks(0x0),
 175   fptTracksCut(0x0),
 176   fetaTracksCut(0x0),
 177   fphiTracksCut(0x0),
 178   fdetaTracksCut(0x0),
 179   fdphiTracksCut(0x0),
 180   fetaphiptTracksCut(0x0),
 181   fetaphiTracksCut(0x0),
 182   fdetadphiTracksCut(0x0),
 183   fNPtTrig(0x0),
 184   fNPtTrigCut(0x0),
 185   fvertexXY(0x0),
 186   fvertexZ(0x0),
 187   fEvtMult(0x0),
 188   fEvtMultvsJetPt(0x0),
 189   fPtvsEtaJet(0x0),
 190   fNpvsEtaJet(0x0),
 191   fNpevtvsEtaJet(0x0),
 192   fPtvsPtJet(0x0),
 193   fNpvsPtJet(0x0),
 194   fNpevtvsPtJet(0x0),
 195   fPtvsPtJet1D(0x0),
 196   fNpvsPtJet1D(0x0),
 197   fNpevtvsPtJet1D(0x0),
 198   fptLeadingJet(0x0),
 199   fetaLeadingJet(0x0),
 200   fphiLeadingJet(0x0),
 201   fptJet(0x0),
 202   fetaJet(0x0),
 203   fphiJet(0x0),
 204   fHistList(0x0),
 205   fNBadRuns(0),
 206   fNBadRunsH(0x0)
 207  {
 208   //
 209   // Default constructor
 210   //
 211   /*
 212   for(int i = 0;i < kMaxStep*2;++i){
 213     fhnJetContainer[i] = 0;
 214   }
 215   */
 216 //   for(int ij  = 0;ij<kMaxJets;++ij){
 217 //     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 218 //     fh1Eta[ij] = fh1Phi[ij] = 0;
 219 //   }
 220 }
 221
 222
 223 //#######################################################################
 224 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const char* name):
 225   AliAnalysisTaskSE(name),
 226   fJetHeaderRec(0x0),
 227   fJetHeaderGen(0x0),
 228   fAOD(0x0),
 229   fBranchRec(""),
 230   fBranchGen(""),
 231   fUseAODInput(kFALSE),
 232   fUseAODJetInput(kFALSE),
 233   fUseAODTrackInput(kFALSE),
 234   fUseAODMCInput(kFALSE),
 235   fUseGlobalSelection(kFALSE),
 236   fUseExternalWeightOnly(kFALSE),
 237   fLimitGenJetEta(kFALSE),
 238   fFilterMask(0),
 239   fAnalysisType(0),
 240   fTrackTypeRec(kTrackUndef),
 241   fTrackTypeGen(kTrackUndef),
 242   fAvgTrials(1.),
 243   fExternalWeight(1.),
 244   fRecEtaWindow(0.5),
 245   fR(0.),
 246   fdRdNdxi(0.7),
 247   fPartPtCut(0.),
 248   fEfactor(1.),
 249   fNff(5),
 250   fNim(5),
 251   fList(0x0),
 252   fGlobVar(1),
 253   fCDFCut(1),
 254   // Number of energy bins
 255   fnEBin(6),
 256   fEmin(10.),
 257   fEmax(70.),
 258   fnEInterval(6),
 259   // Number of radius bins
 260   fnRBin(10),
 261   fRmin(0.1),
 262   fRmax(1.),
 263   fnRInterval(9),
 264   // eta histograms
 265   fnEtaHBin(50),
 266   fEtaHBinMin(-0.9),
 267   fEtaHBinMax(0.9),
 268   // phi histograms
 269   fnPhiHBin(60),
 270   fPhiHBinMin(0.),
 271   fPhiHBinMax(2*TMath::Pi()),
 272   // pt histograms
 273   fnPtHBin(300),
 274   fPtHBinMin(0.),
 275   fPtHBinMax(300.),
 276   // E histograms
 277   fnEHBin(300),
 278   fEHBinMin(0.),
 279   fEHBinMax(300.),
 280   // Xi histograms
 281   fnXiHBin(27),
 282   fXiHBinMin(0.),
 283   fXiHBinMax(9.),
 284   // Pthad histograms
 285   fnPthadHBin(60),
 286   fPthadHBinMin(0.),
 287   fPthadHBinMax(30.),
 288   // z histograms
 289   fnZHBin(30),
 290   fZHBinMin(0.),
 291   fZHBinMax(1.),
 292   fnPtTrigBin(10),
 293   // theta histograms
 294   fnThetaHBin(200),
 295   fThetaHBinMin(-0.5),
 296   fThetaHBinMax(1.5),
 297   fnCosThetaHBin(100),
 298   fcosThetaHBinMin(0.),
 299   fcosThetaHBinMax(1.),
 300   // kT histograms
 301   fnkTHBin(25),
 302   fkTHBinMin(0.),
 303   fkTHBinMax(5.),
 304   // kT histograms
 305   fnRHBin(10),
 306   fRHBinMin(0.),
 307   fRHBinMax(1.),
 308   //Histograms
 309   fEtaMonoJet1H(0x0),
 310   fPhiMonoJet1H(0x0),
 311   fPtMonoJet1H(0x0),
 312   fEMonoJet1H(0x0),
 313   fdNdXiMonoJet1H(0x0),
 314   fdNdPtMonoJet1H(0x0),
 315   fdNdZMonoJet1H(0x0),
 316   fdNdThetaMonoJet1H(0x0),
 317   fdNdcosThetaMonoJet1H(0x0),
 318   fdNdkTMonoJet1H(0x0),
 319   fdNdpTvsZMonoJet1H(0x0),
 320   fShapeMonoJet1H(0x0),
 321   fNMonoJet1sH(0x0),
 322   fThetaPtPartMonoJet1H(0x0),
 323   fcosThetaPtPartMonoJet1H(0x0),
 324   fkTPtPartMonoJet1H(0x0),
 325   fThetaPtJetMonoJet1H(0x0),
 326   fcosThetaPtJetMonoJet1H(0x0),
 327   fkTPtJetMonoJet1H(0x0),
 328   fpTPtJetMonoJet1H(0x0),
 329   farrayEmin(0x0),
 330   farrayEmax(0x0),
 331   farrayRadii(0x0),
 332   farrayPtTrigmin(0x0),
 333   farrayPtTrigmax(0x0),
 334   // Track control plots
 335   fptAllTracks(0x0),
 336   fetaAllTracks(0x0),
 337   fphiAllTracks(0x0),
 338   fetaphiptAllTracks(0x0),
 339   fetaphiAllTracks(0x0),
 340   fptAllTracksCut(0x0),
 341   fetaAllTracksCut(0x0),
 342   fphiAllTracksCut(0x0),
 343   fetaphiptAllTracksCut(0x0),
 344   fetaphiAllTracksCut(0x0),
 345   fptTracks(0x0),
 346   fetaTracks(0x0),
 347   fphiTracks(0x0),
 348   fdetaTracks(0x0),
 349   fdphiTracks(0x0),
 350   fetaphiptTracks(0x0),
 351   fetaphiTracks(0x0),
 352   fdetadphiTracks(0x0),
 353   fptTracksCut(0x0),
 354   fetaTracksCut(0x0),
 355   fphiTracksCut(0x0),
 356   fdetaTracksCut(0x0),
 357   fdphiTracksCut(0x0),
 358   fetaphiptTracksCut(0x0),
 359   fetaphiTracksCut(0x0),
 360   fdetadphiTracksCut(0x0),
 361   fNPtTrig(0x0),
 362   fNPtTrigCut(0x0),
 363   fvertexXY(0x0),
 364   fvertexZ(0x0),
 365   fEvtMult(0x0),
 366   fEvtMultvsJetPt(0x0),
 367   fPtvsEtaJet(0x0),
 368   fNpvsEtaJet(0x0),
 369   fNpevtvsEtaJet(0x0),
 370   fPtvsPtJet(0x0),
 371   fNpvsPtJet(0x0),
 372   fNpevtvsPtJet(0x0),
 373   fPtvsPtJet1D(0x0),
 374   fNpvsPtJet1D(0x0),
 375   fNpevtvsPtJet1D(0x0),
 376   fptLeadingJet(0x0),
 377   fetaLeadingJet(0x0),
 378   fphiLeadingJet(0x0),
 379   fptJet(0x0),
 380   fetaJet(0x0),
 381   fphiJet(0x0),
 382   fHistList(0x0),
 383   fNBadRuns(0),
 384   fNBadRunsH(0x0)
 385 {
 386   //
 387   // Default constructor
 388   //
 389   /*
 390   for(int i = 0;i < kMaxStep*2;++i){
 391     fhnJetContainer[i] = 0;
 392   }
 393   */
 394 //   for(int ij  = 0;ij<kMaxJets;++ij){
 395 //     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 396 //     fh1Eta[ij] = fh1Phi[ij] = 0;
 397 //   }
 398
 399   DefineOutput(1, TList::Class());
 400 }
 401
 402 //////////////////////////////////////////////////////////////////////////////
 403
 404 Bool_t AliAnalysisTaskFragmentationFunction::Notify()
 405 {
 406   //
 407   // Implemented Notify() to read the cross sections
 408   // and number of trials from pyxsec.root
 409   //
 410
 411 //   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 412 //   UInt_t ntrials  = 0;
 413 //   Float_t ftrials  = 0;
 414 //   if(tree){
 415 //     TFile *curfile = tree->GetCurrentFile();
 416 //     if (!curfile) {
 417 //       Error("Notify","No current file");
 418 //       return kFALSE;
 419 //     }
 420
 421 //     if(!fh1Xsec||!fh1Trials){
 422 //       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 423 //       return kFALSE;
 424 //     }
 425
 426 //     TString fileName(curfile->GetName());
 427 //     if(fileName.Contains("AliESDs.root")){
 428 //         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
 429 //     }
 430 //     else if(fileName.Contains("AliAOD.root")){
 431 //         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
 432 //     }
 433 //     else if(fileName.Contains("AliAODs.root")){
 434 //         fileName.ReplaceAll("AliAODs.root", "");
 435 //     }
 436 //     else if(fileName.Contains("galice.root")){
 437 //         // for running with galice and kinematics alone...
 438 //         fileName.ReplaceAll("galice.root", "pyxsec.root");
 439 //     }
 440 //     TFile *fxsec = TFile::Open(fileName.Data());
 441 //     if(!fxsec){
 442 //       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
 443 //       // no a severe condition
 444 //       return kTRUE;
 445 //     }
 446 //     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 447 //     if(!xtree){
 448 //       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
 449 //     }
 450 //     xtree->SetBranchAddress("xsection",&fXsection);
 451 //     xtree->SetBranchAddress("ntrials",&ntrials);
 452 //     ftrials = ntrials;
 453 //     xtree->GetEntry(0);
 454
 455 //     fh1Xsec->Fill("<#sigma>",fXsection);
 456 //     fh1Trials->Fill("#sum{ntrials}",ftrials);
 457
 458 //   }
 459
 460   return kTRUE;
 461
 462 }
 463
 464 //////////////////////////////////////////////////////////////////////////////
 465 //////////////////////////////////////////////////////////////////////////////
 466
 467 void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
 468 {
 469   //
 470   // Create the output container
 471   //
 472
 473   //**** Connect the AOD
 474   if(fUseAODInput) // Use AODs as input not ESDs
 475   {
 476     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 477     if(!fAOD)
 478     {
 479       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 480       return;
 481     }
 482
 483     // fetch the header
 484     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 485     if(!fJetHeaderRec)
 486     {
 487       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
 488     }
 489   }
 490
 491
 492   else // Use the AOD on the flight
 493   {
 494     //  assume that the AOD is in the general output...
 495     fAOD  = AODEvent();
 496     if(!fAOD)
 497     {
 498       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 499       return;
 500     }
 501
 502     //    ((TList*)OutputTree()->GetUserInfo())->Dump();
 503     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 504     if(!fJetHeaderRec)
 505     {
 506       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
 507     }
 508     else
 509     {
 510       if(fDebug>10)fJetHeaderRec->Dump();
 511     }
 512   }
 513
 514 ////////////////////////////
 515
 516   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 517
 518   OpenFile(1);
 519   if(!fHistList)fHistList = new TList();
 520   fHistList->SetOwner(kTRUE);
 521   Bool_t oldStatus = TH1::AddDirectoryStatus();
 522   TH1::AddDirectory(kFALSE);
 523
 524 //////////////////////////////////////////////////
 525 //////////// HISTOGRAM DECLARATION ///////////////
 526 //////////////////////////////////////////////////
 527
 528   DefineJetH();
 529
 530 //////////////////////////////////////////////////
 531 ////////////// HISTOGRAM SAVING //////////////////
 532 //////////////////////////////////////////////////
 533
 534   for (Int_t i3 = 0; i3 < fnEBin; i3++)
 535   {
 536     fHistList->Add(fEtaMonoJet1H[i3]);
 537     fHistList->Add(fPhiMonoJet1H[i3]);
 538     fHistList->Add(fPtMonoJet1H[i3]);
 539     fHistList->Add(fEMonoJet1H[i3]);
 540
 541     for(Int_t i4 = 0; i4 < fnRBin; i4++)
 542     {
 543       fHistList->Add(fdNdXiMonoJet1H[i3][i4]);
 544       fHistList->Add(fdNdPtMonoJet1H[i3][i4]);
 545       fHistList->Add(fdNdZMonoJet1H[i3][i4]);
 546       fHistList->Add(fNMonoJet1sH[i3][i4]);
 547     }
 548   }
 549
 550   // Theta, kT particles/jet
 551   for (Int_t i3 = 0; i3 < fnEBin; i3++)
 552   {
 553     for(Int_t i4 = 0; i4 < fnRBin; i4++)
 554     {
 555       fHistList->Add(fdNdThetaMonoJet1H[i3][i4]);
 556       fHistList->Add(fdNdcosThetaMonoJet1H[i3][i4]);
 557       fHistList->Add(fdNdkTMonoJet1H[i3][i4]);
 558       fHistList->Add(fdNdpTvsZMonoJet1H[i3][i4]);
 559       fHistList->Add(fShapeMonoJet1H[i3][i4]);
 560
 561       fHistList->Add(fThetaPtPartMonoJet1H[i3][i4]);
 562       fHistList->Add(fcosThetaPtPartMonoJet1H[i3][i4]);
 563       fHistList->Add(fkTPtPartMonoJet1H[i3][i4]);
 564       fHistList->Add(fThetaPtJetMonoJet1H[i3][i4]);
 565       fHistList->Add(fcosThetaPtJetMonoJet1H[i3][i4]);
 566       fHistList->Add(fkTPtJetMonoJet1H[i3][i4]);
 567       fHistList->Add(fpTPtJetMonoJet1H[i3][i4]);
 568     }
 569   }
 570
 571   // Track QA - Correlations
 572   for (Int_t iPtBin=0; iPtBin<fnPtTrigBin; iPtBin++)
 573     {
 574       fHistList->Add(fptTracks[iPtBin]);
 575       fHistList->Add(fetaTracks[iPtBin]);
 576       fHistList->Add(fphiTracks[iPtBin]);
 577       fHistList->Add(fdetaTracks[iPtBin]);
 578       fHistList->Add(fdphiTracks[iPtBin]);
 579       fHistList->Add(fetaphiptTracks[iPtBin]);
 580       fHistList->Add(fetaphiTracks[iPtBin]);
 581       fHistList->Add(fdetadphiTracks[iPtBin]);
 582       fHistList->Add(fptTracksCut[iPtBin]);
 583       fHistList->Add(fetaTracksCut[iPtBin]);
 584       fHistList->Add(fphiTracksCut[iPtBin]);
 585       fHistList->Add(fdetaTracksCut[iPtBin]);
 586       fHistList->Add(fdphiTracksCut[iPtBin]);
 587       fHistList->Add(fetaphiptTracksCut[iPtBin]);
 588       fHistList->Add(fetaphiTracksCut[iPtBin]);
 589       fHistList->Add(fdetadphiTracksCut[iPtBin]);
 590       fHistList->Add(fNPtTrig[iPtBin]);
 591       fHistList->Add(fNPtTrigCut[iPtBin]);
 592     }
 593
 594   // Track QA
 595   fHistList->Add(fptAllTracks);
 596   fHistList->Add(fetaAllTracks);
 597   fHistList->Add(fphiAllTracks);
 598   fHistList->Add(fetaphiptAllTracks);
 599   fHistList->Add(fetaphiAllTracks);
 600   fHistList->Add(fptAllTracksCut);
 601   fHistList->Add(fetaAllTracksCut);
 602   fHistList->Add(fphiAllTracksCut);
 603   fHistList->Add(fetaphiptAllTracksCut);
 604   fHistList->Add(fetaphiAllTracksCut);
 605
 606   // Event caracterisation QA
 607   fHistList->Add(fvertexXY);
 608   fHistList->Add(fvertexZ);
 609   fHistList->Add(fEvtMult);
 610   fHistList->Add(fEvtMultvsJetPt);
 611   fHistList->Add(fPtvsEtaJet);
 612   fHistList->Add(fNpvsEtaJet);
 613   fHistList->Add(fNpevtvsEtaJet);
 614   fHistList->Add(fPtvsPtJet);
 615   fHistList->Add(fNpvsPtJet);
 616   fHistList->Add(fNpevtvsPtJet);
 617   fHistList->Add(fPtvsPtJet1D);
 618   fHistList->Add(fNpvsPtJet1D);
 619   fHistList->Add(fNpevtvsPtJet1D);
 620   fHistList->Add(fptLeadingJet);
 621   fHistList->Add(fetaLeadingJet);
 622   fHistList->Add(fphiLeadingJet);
 623   fHistList->Add(fptJet);
 624   fHistList->Add(fetaJet);
 625   fHistList->Add(fphiJet);
 626   fHistList->Add(fNBadRunsH);
 627
 628 //////////////////////////////////////////////////
 629 ///////// END OF HISTOGRAM DECLARATION ///////////
 630 //////////////////////////////////////////////////
 631
 632   // =========== Switch on Sumw2 for all histos ===========
 633   for (Int_t i=0; i<fHistList->GetEntries(); ++i)
 634   {
 635     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 636     if (h1)
 637     {
 638       // Printf("%s ",h1->GetName());
 639       h1->Sumw2();
 640       continue;
 641     }
 642   }
 643
 644   TH1::AddDirectory(oldStatus);
 645
 646   /*
 647     if (fDebug > 1) printf("AnalysisTaskDiJets::CreateOutPutData() \n");
 648     fDiJets = new TClonesArray("AliAODDiJet", 0);
 649     fDiJets->SetName("Dijets");
 650     AddAODBranch("TClonesArray", &fDiJets);
 651   */
 652
 653
 654 }
 655
 656 //////////////////////////////////////////////////////////////////////////////
 657 //////////////////////////////////////////////////////////////////////////////
 658
 659 void AliAnalysisTaskFragmentationFunction::Init()
 660 {
 661   //
 662   // Initialization
 663   //
 664
 665   Printf(">>> AnalysisTaskFragmentationFunction::Init() debug level %d\n",fDebug);
 666   if (fDebug > 1) printf("AnalysisTaskDiJets::Init() \n");
 667 }
 668
 669 /////////////////////////////////////////////////////////////////////////////
 670 /////////////////////////////////////////////////////////////////////////////
 671
 672 void AliAnalysisTaskFragmentationFunction::UserExec(Option_t */*option*/)
 673 {
 674   //
 675   // Execute analysis for current event
 676   //
 677
 678   //****
 679   //**** Check of input data
 680   //****
 681
 682   printf("Analysing event # %5d\n", (Int_t) fEntry);
 683   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 684
 685   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 686   if(!aodH)
 687   {
 688     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 689     return;
 690   }
 691
 692   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 693   if(!aodRecJets)
 694   {
 695     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 696     return;
 697   }
 698   if (fDebug > 10) Printf("%s:%d",(char*)__FILE__,__LINE__);
 699
 700   //****
 701   //**** Check of primary vertex
 702   //****
 703   AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
 704   if( (!pvtx) ||
 705       (pvtx->GetZ()<-10. || pvtx->GetZ()>10.) ||
 706       (pvtx->GetNContributors()<0) )
 707     {
 708         fNBadRuns++;
 709         fNBadRunsH->Fill(0.5);
 710         return;
 711     }
 712
 713   //****
 714   //**** Check number of tracks
 715   //****
 716   TClonesArray* tracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 717
 718   //****
 719   //**** Declaration of arrays and variables
 720   //****
 721   // We use static array, not to fragment the memory
 722   AliAODJet recJets[kMaxJets];
 723   Int_t nRecJets       = 0;
 724   Int_t nTracks = 0;
 725
 726 //////////////////////////////////////////////////
 727 ///////// Get the reconstructed jets /////////////
 728 //////////////////////////////////////////////////
 729
 730   nRecJets = aodRecJets->GetEntries();
 731   nRecJets = TMath::Min(nRecJets,kMaxJets);
 732   nTracks  = fAOD->GetNumberOfTracks();
 733
 734   for(int ir = 0;ir < nRecJets;++ir)
 735   {
 736     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 737     if(!tmp)continue;
 738     recJets[ir] = *tmp;
 739     cout << "recJets[" << ir << "].Eta(): " << recJets[ir].Eta() << ", recJets[" << ir <<"].Phi(): " << recJets[ir].Phi() << ", recJets[" << ir << "].E(): " << recJets[ir].E() << endl;
 740   }
 741   if(nRecJets>1) {
 742     Float_t detaJ = recJets[0].Eta() - recJets[1].Eta();
 743     Float_t dphiJ = recJets[0].Phi() - recJets[1].Phi();
 744     Float_t detJ = recJets[0].Pt() - recJets[1].Pt();
 745     cout << "detaJ: " << detaJ << ", dphiJ: " << dphiJ << ", detJ: " << detJ << endl;
 746   }
 747
 748   // Get vertex information
 749   fvertexXY->Fill(pvtx->GetX(),pvtx->GetY());
 750   fvertexZ->Fill(pvtx->GetZ());
 751
 752   //////////////////////////////////////////////////
 753   ////////////// TRACK QUALITY ASSURANCE ///////////           TO BE OPTIMISED!!
 754   //////////////////////////////////////////////////
 755   Int_t evtMult = 0;
 756   for(Int_t it=0; it<nTracks; it++)
 757   {
 758     AliAODTrack* aodTrack = (AliAODTrack*)tracks->At(it);
 759     Float_t etaT = aodTrack->Eta();
 760     Float_t phiT = aodTrack->Phi();
 761     Float_t ptT = aodTrack->Pt();
 762     phiT = ((phiT < 0) ? phiT + 2 * TMath::Pi() : phiT);
 763     //    cout << "etaT: " << etaT << ", phiT: " << phiT << endl;
 764     fptAllTracks->Fill(ptT);
 765     fetaAllTracks->Fill(etaT);
 766     fphiAllTracks->Fill(phiT);
 767     fetaphiptAllTracks->Fill(etaT,phiT,ptT);
 768     fetaphiAllTracks->Fill(etaT,phiT,1);
 769     UInt_t status = aodTrack->GetStatus();
 770
 771     for(Int_t i=0; i<fnPtTrigBin; i++)
 772     {
 773       if(ptT>=farrayPtTrigmin[i] && ptT<farrayPtTrigmax[i])
 774       {
 775         fptTracks[i]->Fill(ptT);
 776         fetaTracks[i]->Fill(etaT);
 777         fphiTracks[i]->Fill(phiT);
 778         fNPtTrig[i]->Fill(0.5);
 779         // Compute deta/dphi
 780         Float_t etaT2 = 0.;
 781         Float_t phiT2 = 0.;
 782         Float_t ptT2 = 0.;
 783         Float_t detaT = 0.;
 784         Float_t dphiT = 0.;
 785         for(Int_t it2 = 0; it2< nTracks; it2++)
 786         {
 787            AliAODTrack* aodTrack2 = (AliAODTrack*)tracks->At(it2);
 788            etaT2 = aodTrack2->Eta(); phiT2 = aodTrack2->Phi(); ptT2 = aodTrack2->Pt();
 789            phiT2 = ((phiT2 < 0) ? phiT2 + 2 * TMath::Pi() : phiT2);
 790            //         cout << "etaT2: " << etaT2 << ", phiT2: " << phiT2 << endl;
 791            if(ptT2 > 2.*i+4.) continue;
 792            if(it2==it) continue;
 793            detaT = etaT - etaT2;
 794            dphiT = phiT - phiT2;
 795            if (dphiT  >   TMath::Pi()) dphiT = (-TMath::Pi() +TMath::Abs(dphiT - TMath::Pi()));
 796            if (dphiT  < -1.0*TMath::Pi())  dphiT = (TMath::Pi() -  TMath::Abs(dphiT + TMath::Pi()));
 797
 798            fdetaTracks[i]->Fill(detaT);
 799            fdphiTracks[i]->Fill(dphiT);
 800            fdetadphiTracks[i]->Fill(detaT,dphiT,1);
 801         }
 802         fetaphiptTracks[i]->Fill(etaT,phiT,ptT);
 803         fetaphiTracks[i]->Fill(etaT,phiT,1);
 804       }
 805     } // End loop over trigger ranges
 806
 807     if (status == 0) continue;
 808     if((fFilterMask>0)&&!(aodTrack->TestFilterBit(fFilterMask))) continue;
 809     fptAllTracksCut->Fill(ptT);
 810     fetaAllTracksCut->Fill(etaT);
 811     fphiAllTracksCut->Fill(phiT);
 812     fetaphiptAllTracksCut->Fill(etaT,phiT,ptT);
 813     fetaphiAllTracksCut->Fill(etaT,phiT,1);
 814     if(ptT > 0.150 && TMath::Abs(etaT) < 0.9) evtMult++;
 815   } // end loop over tracks
 816   fEvtMult->Fill(evtMult);
 817
 818 //////////////////////////////////////////////////
 819 ///////////////// MONOJET PART ///////////////////
 820 //////////////////////////////////////////////////
 821
 822   if (nRecJets == 0) return;
 823
 824   Double_t jetEnergy = recJets[0].E();
 825   Int_t    goodBin   = 999;
 826
 827   for (Int_t i1 = 0; i1 < fnEBin; i1++)
 828   {
 829     if (jetEnergy < farrayEmax[i1] && jetEnergy >= farrayEmin[i1])
 830     {
 831       goodBin = i1;
 832       continue;
 833     }
 834   }
 835
 836   fptLeadingJet->Fill(recJets[0].Pt());
 837   fetaLeadingJet->Fill(recJets[0].Eta());
 838   fphiLeadingJet->Fill(recJets[0].Phi());
 839
 840   for(Int_t ij=0; ij<nRecJets; ij++)
 841   {
 842     fptJet->Fill(recJets[ij].Pt());
 843     fetaJet->Fill(recJets[ij].Eta());
 844     fphiJet->Fill(recJets[ij].Phi());
 845   }
 846
 847   // Get track ref
 848   TRefArray* ref = recJets[0].GetRefTracks();
 849   for(Int_t it=0; it<ref->GetEntries(); it++)
 850   {
 851     Float_t ptTrack = ((AliVTrack*)ref->At(it))->Pt();
 852     fPtvsEtaJet->Fill(recJets[0].Eta(),ptTrack);
 853     fNpvsEtaJet->Fill(recJets[0].Eta(),ref->GetEntries());
 854     fNpevtvsEtaJet->Fill(recJets[0].Eta(),evtMult);
 855     fPtvsPtJet->Fill(recJets[0].Pt(),ptTrack);
 856     fNpvsPtJet->Fill(recJets[0].Pt(),ref->GetEntries());
 857     fNpevtvsPtJet->Fill(recJets[0].Pt(),evtMult);
 858     fPtvsPtJet1D->Fill(recJets[0].Pt(),ptTrack);
 859     fNpvsPtJet1D->Fill(recJets[0].Pt(),ref->GetEntries());
 860     fNpevtvsPtJet1D->Fill(recJets[0].Pt(),evtMult);
 861   }
 862
 863   FillMonoJetH(goodBin, recJets, tracks);
 864
 865   //////////////////////////////////////////////////
 866   ////////////////// DIJET PART ////////////////////
 867   //////////////////////////////////////////////////
 868
 869   // UNDER CONSTRUCTION
 870
 871   PostData(1, fHistList);
 872 }
 873
 874 //#######################################################################
 875 void AliAnalysisTaskFragmentationFunction::Terminate(Option_t */*option*/)
 876 {
 877 // Terminate analysis
 878 //
 879     if (fDebug > 1) printf("AnalysisDiJets: Terminate() \n");
 880     printf("Number of events with vertex out of bound: %d", fNBadRuns);
 881 }
 882
 883 //////////////////////////////////////////////////////////////////////////////
 884 //////////////////////////////////////////////////////////////////////////////
 885
 886 void AliAnalysisTaskFragmentationFunction::DefineJetH()
 887 {
 888
 889   /////////////////////////////////////// HISTOGRAMS FIRST JET
 890   fEtaMonoJet1H   = new TH1F*[fnEBin+1];
 891   fPhiMonoJet1H   = new TH1F*[fnEBin+1];
 892   fPtMonoJet1H    = new TH1F*[fnEBin+1];
 893   fEMonoJet1H     = new TH1F*[fnEBin+1];
 894
 895   fdNdXiMonoJet1H = new TH1F**[fnEBin+1];
 896   fdNdPtMonoJet1H = new TH1F**[fnEBin+1];
 897   fdNdZMonoJet1H  = new TH1F**[fnEBin+1];
 898   fdNdThetaMonoJet1H    = new TH1F**[fnEBin+1];
 899   fdNdcosThetaMonoJet1H    = new TH1F**[fnEBin+1];
 900   fdNdkTMonoJet1H       = new TH1F**[fnEBin+1];
 901   fdNdpTvsZMonoJet1H    = new TH1F**[fnEBin+1];
 902   fShapeMonoJet1H       = new TH1F**[fnEBin+1];
 903   fNMonoJet1sH          = new TH1F**[fnEBin+1];
 904
 905   fThetaPtPartMonoJet1H = new TH2F**[fnEBin+1];
 906   fcosThetaPtPartMonoJet1H = new TH2F**[fnEBin+1];
 907   fkTPtPartMonoJet1H    = new TH2F**[fnEBin+1];
 908   fThetaPtJetMonoJet1H = new TH2F**[fnEBin+1];
 909   fcosThetaPtJetMonoJet1H = new TH2F**[fnEBin+1];
 910   fkTPtJetMonoJet1H    = new TH2F**[fnEBin+1];
 911   fpTPtJetMonoJet1H    = new TH2F**[fnEBin+1];
 912
 913   for (Int_t iEbin=0;iEbin<fnEBin+1;iEbin++)
 914   {
 915     fdNdXiMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
 916     fdNdPtMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
 917     fdNdZMonoJet1H[iEbin]        = new TH1F*[fnRBin+1];
 918     fdNdThetaMonoJet1H[iEbin]    = new TH1F*[fnRBin+1];
 919     fdNdcosThetaMonoJet1H[iEbin] = new TH1F*[fnRBin+1];
 920     fdNdkTMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
 921     fdNdpTvsZMonoJet1H[iEbin]    = new TH1F*[fnRBin+1];
 922     fShapeMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
 923     fNMonoJet1sH[iEbin]          = new TH1F*[fnRBin+1];
 924
 925     fThetaPtPartMonoJet1H[iEbin]    = new TH2F*[fnRBin+1];
 926     fcosThetaPtPartMonoJet1H[iEbin] = new TH2F*[fnRBin+1];
 927     fkTPtPartMonoJet1H[iEbin]       = new TH2F*[fnRBin+1];
 928     fThetaPtJetMonoJet1H[iEbin]     = new TH2F*[fnRBin+1];
 929     fcosThetaPtJetMonoJet1H[iEbin]  = new TH2F*[fnRBin+1];
 930     fkTPtJetMonoJet1H[iEbin]        = new TH2F*[fnRBin+1];
 931     fpTPtJetMonoJet1H[iEbin]        = new TH2F*[fnRBin+1];
 932   }
 933
 934   for (Int_t iEbin=0;iEbin<fnEBin+1;iEbin++)
 935   {
 936     fEtaMonoJet1H[iEbin]    = 0;
 937     fPhiMonoJet1H[iEbin]    = 0;
 938     fPtMonoJet1H[iEbin]     = 0;
 939     fEMonoJet1H[iEbin]      = 0;
 940
 941     for (Int_t iRbin=0;iRbin<fnRBin+1;iRbin++)
 942     {
 943       fdNdXiMonoJet1H[iEbin][iRbin]       = 0;
 944       fdNdPtMonoJet1H[iEbin][iRbin]       = 0;
 945       fdNdZMonoJet1H[iEbin][iRbin]        = 0;
 946       fdNdThetaMonoJet1H[iEbin][iRbin]    = 0;
 947       fdNdcosThetaMonoJet1H[iEbin][iRbin] = 0;
 948       fdNdkTMonoJet1H[iEbin][iRbin]       = 0;
 949       fdNdpTvsZMonoJet1H[iEbin][iRbin]    = 0;
 950       fShapeMonoJet1H[iEbin][iRbin]       = 0;
 951       fNMonoJet1sH[iEbin][iRbin]          = 0;
 952
 953       fThetaPtPartMonoJet1H[iEbin][iRbin]    = 0;
 954       fcosThetaPtPartMonoJet1H[iEbin][iRbin] = 0;
 955       fkTPtPartMonoJet1H[iEbin][iRbin]       = 0;
 956       fThetaPtJetMonoJet1H[iEbin][iRbin]     = 0;
 957       fcosThetaPtJetMonoJet1H[iEbin][iRbin]  = 0;
 958       fkTPtJetMonoJet1H[iEbin][iRbin]        = 0;
 959       fpTPtJetMonoJet1H[iEbin][iRbin]        = 0;
 960     }
 961   }
 962
 963   fptTracks          = new TH1F*[fnPtTrigBin+1];
 964   fetaTracks         = new TH1F*[fnPtTrigBin+1];
 965   fphiTracks         = new TH1F*[fnPtTrigBin+1];
 966   fdetaTracks        = new TH1F*[fnPtTrigBin+1];
 967   fdphiTracks        = new TH1F*[fnPtTrigBin+1];
 968   fetaphiptTracks    = new TH2F*[fnPtTrigBin+1];
 969   fetaphiTracks      = new TH2F*[fnPtTrigBin+1];
 970   fdetadphiTracks    = new TH2F*[fnPtTrigBin+1];
 971   fptTracksCut       = new TH1F*[fnPtTrigBin+1];
 972   fetaTracksCut      = new TH1F*[fnPtTrigBin+1];
 973   fphiTracksCut      = new TH1F*[fnPtTrigBin+1];
 974   fdetaTracksCut     = new TH1F*[fnPtTrigBin+1];
 975   fdphiTracksCut     = new TH1F*[fnPtTrigBin+1];
 976   fetaphiptTracksCut = new TH2F*[fnPtTrigBin+1];
 977   fetaphiTracksCut   = new TH2F*[fnPtTrigBin+1];
 978   fdetadphiTracksCut = new TH2F*[fnPtTrigBin+1];
 979   fNPtTrig           = new TH1F*[fnPtTrigBin+1];
 980   fNPtTrigCut        = new TH1F*[fnPtTrigBin+1];
 981
 982   for(Int_t iPtTrigBin=0; iPtTrigBin<fnPtTrigBin; iPtTrigBin++)
 983   {
 984     fptTracks[iPtTrigBin]          = 0;
 985     fetaTracks[iPtTrigBin]         = 0;
 986     fphiTracks[iPtTrigBin]         = 0;
 987     fdetaTracks[iPtTrigBin]        = 0;
 988     fdphiTracks[iPtTrigBin]        = 0;
 989     fetaphiptTracks[iPtTrigBin]    = 0;
 990     fetaphiTracks[iPtTrigBin]      = 0;
 991     fdetadphiTracks[iPtTrigBin]    = 0;
 992     fptTracksCut[iPtTrigBin]       = 0;
 993     fetaTracksCut[iPtTrigBin]      = 0;
 994     fphiTracksCut[iPtTrigBin]      = 0;
 995     fdetaTracksCut[iPtTrigBin]     = 0;
 996     fdphiTracksCut[iPtTrigBin]     = 0;
 997     fetaphiptTracksCut[iPtTrigBin] = 0;
 998     fetaphiTracksCut[iPtTrigBin]   = 0;
 999     fdetadphiTracksCut[iPtTrigBin] = 0;
1000     fNPtTrig[iPtTrigBin]           = 0;
1001     fNPtTrigCut[iPtTrigBin]        = 0;
1002   }
1003
1004   farrayEmin = new Double_t[fnEBin];
1005   farrayEmax = new Double_t[fnEBin];
1006
1007   farrayPtTrigmin = new Double_t[fnPtTrigBin];
1008   farrayPtTrigmax = new Double_t[fnPtTrigBin];
1009
1010   farrayRadii = new Double_t[fnRBin];
1011
1012   Double_t Emin    = 0;
1013   Double_t Emax    = 0;
1014   Double_t pasE    = (Double_t)((fEmax-fEmin)/fnEInterval);
1015   TString energy;
1016   TString energy2;
1017
1018   Double_t R    = 0.;
1019   Double_t pasR = (Double_t)((fRmax-fRmin)/fnRInterval);
1020   TString radius;
1021
1022   for (Int_t i = 0; i < fnEBin; i++)
1023   {
1024     Emin       = 0;
1025     Emax       = 0;
1026
1027     if (i==0) Emin = fEmin;
1028     if (i!=0) Emin = fEmin + pasE*i;
1029     Emax = Emin+pasE;
1030     energy2 = "E_{jet1} : ";
1031     energy = "E_{jet} : ";
1032     energy += Emin;
1033     energy2 += Emin;
1034     energy += "-";
1035     energy2 += "-";
1036     energy +=Emax;
1037     energy2 += Emax;
1038     energy += "GeV";
1039     energy2 += "GeV";
1040
1041     farrayEmin[i]      = Emin;
1042     farrayEmax[i]      = Emax;
1043
1044     for (Int_t j = 0; j < fnRBin; j++)
1045     {
1046
1047       if (j==0) R = fRmin;
1048       if (j!=0) R = fRmin + pasR*j;
1049       radius = ", R = ";
1050       radius += R;
1051
1052       farrayRadii[j] = R;
1053
1054       fEtaMonoJet1H[i]   = new TH1F("fEtaMonoJet1H,"+energy, "#eta_{jet1},"+energy, fnEtaHBin, fEtaHBinMin, fEtaHBinMax);
1055       fPhiMonoJet1H[i]   = new TH1F("fPhiMonoJet1H,"+energy, "#phi_{jet1},"+energy, fnPhiHBin, fPhiHBinMin, fPhiHBinMax);
1056       fPtMonoJet1H[i]    = new TH1F("fPtMonoJet1H,"+energy, "pT_{jet1},"+energy, fnPtHBin, fPtHBinMin, fPtHBinMax);
1057       fEMonoJet1H[i]     = new TH1F("fEMonoJet1H,"+energy, "E_{jet1},"+energy, fnEHBin, fEHBinMin, fEHBinMax);
1058
1059       fdNdXiMonoJet1H[i][j] = new TH1F("fdNdXiMonoJet1H,"+energy+radius, "dN_{ch}/d#xi,"+energy+radius, fnXiHBin, fXiHBinMin, fXiHBinMax);
1060       fdNdPtMonoJet1H[i][j] = new TH1F("fdNdPtMonoJet1H,"+energy+radius, "dN_{ch}/dPt_{had},"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax);
1061       fdNdZMonoJet1H[i][j]  = new TH1F("fdNdZMonoJet1H,"+energy+radius, "dN_{ch}/dz,"+energy+radius, fnZHBin, fZHBinMin, fZHBinMax);
1062
1063       fdNdThetaMonoJet1H[i][j]    = new TH1F("fdNdThetaMonoJet1H,"+energy+radius, "dN_{ch}/d#Theta,"+energy+radius, fnThetaHBin, fThetaHBinMin, fThetaHBinMax);
1064       fdNdcosThetaMonoJet1H[i][j] = new TH1F("fdNdcosThetaMonoJet1H,"+energy+radius, "dN_{ch}/dcos(#Theta),"+energy+radius, fnCosThetaHBin, fcosThetaHBinMin, fcosThetaHBinMax);
1065       fdNdkTMonoJet1H[i][j]   = new TH1F("fdNdkTMonoJet1H,"+energy+radius, "dN_{ch}/dk_{T},"+energy+radius, fnkTHBin, fkTHBinMin, fkTHBinMax);
1066       fdNdpTvsZMonoJet1H[i][j]= new TH1F("fdNdpTvsZMonoJet1H,"+energy+radius, "dN_{ch}/dk_{T} vs Z,"+energy+radius, fnZHBin, fZHBinMin, fZHBinMax);
1067       fShapeMonoJet1H[i][j]   = new TH1F("fShapeMonoJet1H,"+energy+radius, "E(R=x)/E(R=1)"+energy+radius, fnRHBin, fRHBinMin, fRHBinMax);
1068
1069       fThetaPtPartMonoJet1H[i][j] = new TH2F("fThetaPtPartMonoJet1H,"+energy+radius, "#Theta vs Pt particle,"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax, fnThetaHBin, fThetaHBinMin, fThetaHBinMax);
1070       fcosThetaPtPartMonoJet1H[i][j] = new TH2F("fcosThetaPtPartMonoJet1H,"+energy+radius, "cos(#Theta) vs Pt particle,"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax, fnCosThetaHBin, fcosThetaHBinMin, fcosThetaHBinMax);
1071       fkTPtPartMonoJet1H[i][j] = new TH2F("fkTPtPartMonoJet1H,"+energy+radius, "kT vs Pt particle,"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax, fnkTHBin, fkTHBinMin, fkTHBinMax);
1072       fThetaPtJetMonoJet1H[i][j] = new TH2F("fThetaPtJetMonoJet1H,"+energy+radius, "#Theta vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnThetaHBin, fThetaHBinMin, fThetaHBinMax);
1073       fcosThetaPtJetMonoJet1H[i][j] = new TH2F("fcosThetaPtJetMonoJet1H,"+energy+radius, "cos(#Theta) vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnCosThetaHBin, fcosThetaHBinMin, fcosThetaHBinMax);
1074       fkTPtJetMonoJet1H[i][j] = new TH2F("fkTPtJetMonoJet1H,"+energy+radius, "kT vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnkTHBin, fkTHBinMin, fkTHBinMax);
1075       fpTPtJetMonoJet1H[i][j] = new TH2F("fpTPtJetMonoJet1H,"+energy+radius, "pT vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnkTHBin, fkTHBinMin, fkTHBinMax);
1076
1077       fNMonoJet1sH[i][j]    = new TH1F("fNMonoJet1sH,"+energy+radius, "N_{jets1},"+energy+radius, 1, 0., 1.);
1078       fNBadRunsH = new TH1F("fNBadRunsH","Number of events with Z vertex out of range", 1, 0., 1.);
1079
1080       SetProperties(fEtaMonoJet1H[i], "#eta_{jet1}", "Entries");
1081       SetProperties(fPhiMonoJet1H[i], "#phi_{jet1}", "Entries");
1082       SetProperties(fPtMonoJet1H[i], "p_{Tjet1} (GeV/c)", "Entries");
1083       SetProperties(fEMonoJet1H[i], "E_{jet1} (GeV)", "Entries");
1084
1085       SetProperties(fdNdXiMonoJet1H[i][j], "#xi = ln(E_{jet1}/p_{Thad})", "dN_{had}/d#xi");
1086       SetProperties(fdNdPtMonoJet1H[i][j], "p_{Thad} (GeV/c)", "dN_{had}/dp_{Thad}");
1087       SetProperties(fdNdZMonoJet1H[i][j], "z = (p_{Thad}/E_{jet1})", "dN_{had}/dz");
1088       SetProperties(fdNdThetaMonoJet1H[i][j], "#Theta", "dN_{had}/d#Theta");
1089       SetProperties(fdNdcosThetaMonoJet1H[i][j], "cos(#Theta)", "dN_{had}/dcos(#Theta)");
1090       SetProperties(fdNdkTMonoJet1H[i][j], "k_{Thad}", "dN_{had}/dk_{Thad}");
1091       SetProperties(fdNdpTvsZMonoJet1H[i][j], "z = (p_{Thad}/E_{jet1})", "dN_{had}/dp_{T}");
1092       SetProperties(fShapeMonoJet1H[i][j], "R", "#Psi(R)");
1093       SetProperties(fNMonoJet1sH[i][j], "Bin", "N_{jets1}");
1094
1095       SetProperties(fThetaPtPartMonoJet1H[i][j],"p_{Thad} [GeV/c]","#Theta");
1096       SetProperties(fcosThetaPtPartMonoJet1H[i][j],"p_{Thad} [GeV/c]","cos(#Theta)");
1097       SetProperties(fkTPtPartMonoJet1H[i][j],"p_{Thad} [GeV/c]","k_{Thad}");
1098       SetProperties(fThetaPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "#Theta");
1099       SetProperties(fcosThetaPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "cos(#Theta)");
1100       SetProperties(fkTPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "k_{Thad} [GeV/c]");
1101       SetProperties(fpTPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "p_{Thad} [GeV/c]");
1102     }
1103   }
1104
1105   for(Int_t i=0; i<fnPtTrigBin; i++)
1106   {
1107     if(i==0) farrayPtTrigmin[i] = 1.;
1108     else farrayPtTrigmin[i] = i*5.;
1109     farrayPtTrigmax[i] = i*5+5.;
1110
1111     TString ptTrigRange;
1112     ptTrigRange = "; p_{T} trig range: ";
1113     ptTrigRange += farrayPtTrigmin[i];
1114     ptTrigRange += "-";
1115     ptTrigRange += farrayPtTrigmax[i];
1116     ptTrigRange += " [GeV]";
1117
1118     fptTracks[i] = new TH1F("fptTracks"+ptTrigRange, "Track transverse momentum [GeV]"+ptTrigRange,300,0.,150.);
1119     fetaTracks[i] = new TH1F("fetaTracks"+ptTrigRange, "#eta tracks"+ptTrigRange, 36, -0.9, 0.9);
1120     fphiTracks[i] = new TH1F("fphiTracks"+ptTrigRange, "#phi tracks"+ptTrigRange, 60, 0., 2*TMath::Pi());
1121     fdetaTracks[i] = new TH1F("fdetaTracks"+ptTrigRange, "#Delta #eta tracks"+ptTrigRange,80, -2., 2.);
1122     fdphiTracks[i] = new TH1F("fdphiTracks"+ptTrigRange, "#Delta #phi tracks"+ptTrigRange, 120, -TMath::Pi(), TMath::Pi());
1123     fetaphiptTracks[i] = new TH2F("fetaphiptTracks"+ptTrigRange,"#eta/#phi track p_{T} mapping"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1124     fetaphiTracks[i] = new TH2F("fetaphiTracks"+ptTrigRange,"#eta/#phi track mapping"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1125     fdetadphiTracks[i] = new TH2F("fdetadphiTracks"+ptTrigRange,"#Delta #eta/#Delta #phi track mapping"+ptTrigRange,80, -2., 2., 120, -TMath::Pi(), TMath::Pi());
1126     fptTracksCut[i] = new TH1F("fptTracksCut"+ptTrigRange, "Track transverse momentum after cuts [GeV]"+ptTrigRange,300,0.,150.);
1127     fetaTracksCut[i] = new TH1F("fetaTracksCut"+ptTrigRange, "#eta tracks after cuts"+ptTrigRange, 36, -0.9, 0.9);
1128     fphiTracksCut[i] = new TH1F("fphiTracksCuts"+ptTrigRange, "#phi tracks after cuts"+ptTrigRange, 60, 0., 2*TMath::Pi());
1129     fdetaTracksCut[i] = new TH1F("fdetaTracksCuts"+ptTrigRange, "#Delta #eta tracks after cuts"+ptTrigRange,80, -2., 2.);
1130     fdphiTracksCut[i] = new TH1F("fdphiTracksCuts"+ptTrigRange, "#Delta #phi tracks after cuts"+ptTrigRange, 120, -TMath::Pi(), TMath::Pi());
1131     fetaphiptTracksCut[i] = new TH2F("fetaphiptTracksCuts"+ptTrigRange,"#eta/#phi track p_{T} mapping after cuts"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1132     fetaphiTracksCut[i] = new TH2F("fetaphiTracksCuts"+ptTrigRange,"#eta/#phi track mapping after cuts"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1133     fdetadphiTracksCut[i] = new TH2F("fdetadphiTracksCuts"+ptTrigRange,"#Delta #eta/#Delta #phi track mapping after cuts"+ptTrigRange,80, -2., 2., 120, -TMath::Pi(), TMath::Pi());
1134     fNPtTrig[i] = new TH1F("fNPtTrig"+ptTrigRange,"Number of triggers"+ptTrigRange,1,0.,1.);
1135     fNPtTrigCut[i] = new TH1F("fNPtTrigCut"+ptTrigRange,"Number of triggers after cut"+ptTrigRange,1,0.,1.);
1136
1137     SetProperties(fptTracks[i], "Track p_{T} [GeV]", "dN/dp_{T}");
1138     SetProperties(fetaTracks[i], "Track #eta", "dN/d#eta");
1139     SetProperties(fphiTracks[i], "Track #phi", "dN/d#phi");
1140     SetProperties(fdetaTracks[i], "#eta_{track} - #eta_{trig}", "dN/d#Delta#eta");
1141     SetProperties(fdphiTracks[i], "#phi_{track} - #phi_{trig}", "dN/d#Delta#phi");
1142     SetProperties(fetaphiptTracks[i], "#eta_{track}", "#phi_{track}");
1143     SetProperties(fetaphiTracks[i], "#eta_{track}", "#phi_{track}");
1144     SetProperties(fdetadphiTracks[i], "#Delta #eta_{track}", "#Delta #phi_{track}");
1145     SetProperties(fptTracksCut[i], "p_{T}track [GeV]", "dN/dp_{T}");
1146     SetProperties(fetaTracksCut[i], "#eta_{track}", "dN/d#eta");
1147     SetProperties(fphiTracksCut[i], "#phi_{track}", "dN/d#phi");
1148     SetProperties(fdetaTracksCut[i], "#eta_{track} - #eta_{trig}", "dN/d#Delta#eta");
1149     SetProperties(fdphiTracksCut[i], "#phi_{track} - #phi_{trig}", "dN/d#Delta#phi");
1150     SetProperties(fetaphiptTracksCut[i], "#eta_{track}", "#phi_{track}");
1151     SetProperties(fetaphiTracksCut[i], "#eta_{track}", "#phi_{track}");
1152     SetProperties(fdetadphiTracksCut[i], "#Delta #eta_{track}", "#Delta #phi_{track}");
1153     SetProperties(fNPtTrig[i], "", "Number of triggers");
1154     SetProperties(fNPtTrigCut[i], "", "Number of triggers");
1155
1156   }
1157
1158   fptAllTracks = new TH1F("fptAllTracks", "Track transverse momentum [GeV]",300,0.,150.);
1159   fetaAllTracks = new TH1F("fetaAllTracks", "#eta tracks", 36, -0.9, 0.9);
1160   fphiAllTracks = new TH1F("fphiAllTracks", "#phi tracks", 60, 0., 2*TMath::Pi());
1161   fetaphiptAllTracks = new TH2F("fetaphiptAllTracks","#eta/#phi track p_{T} mapping",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1162   fetaphiAllTracks = new TH2F("fetaphiAllTracks","#eta/#phi track mapping",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1163   fptAllTracksCut = new TH1F("fptAllTracksCut", "Track transverse momentum after cuts [GeV]",300,0.,150.);
1164   fetaAllTracksCut = new TH1F("fetaAllTracksCut", "#eta tracks after cuts", 36, -0.9, 0.9);
1165   fphiAllTracksCut = new TH1F("fphiAllTracksCuts", "#phi tracks after cuts", 60, 0., 2*TMath::Pi());
1166   fetaphiptAllTracksCut = new TH2F("fetaphiptAllTracksCuts","#eta/#phi track p_{T} mapping after cuts",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1167   fetaphiAllTracksCut = new TH2F("fetaphiAllTracksCuts","#eta/#phi track mapping after cuts",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
1168
1169   SetProperties(fptAllTracks, "Track p_{T} [GeV]", "dN/dp_{T}");
1170   SetProperties(fetaAllTracks, "Track #eta", "dN/d#eta");
1171   SetProperties(fphiAllTracks, "Track #phi", "dN/d#phi");
1172   SetProperties(fetaphiptAllTracks, "#eta_{track}", "#phi_{track}");
1173   SetProperties(fetaphiAllTracks, "#eta_{track}", "#phi_{track}");
1174   SetProperties(fptAllTracksCut, "p_{T}track [GeV]", "dN/dp_{T}");
1175   SetProperties(fetaAllTracksCut, "#eta_{track}", "dN/d#eta");
1176   SetProperties(fphiAllTracksCut, "#phi_{track}", "dN/d#phi");
1177   SetProperties(fetaphiptAllTracksCut, "#eta_{track}", "#phi_{track}");
1178   SetProperties(fetaphiAllTracksCut, "#eta_{track}", "#phi_{track}");
1179
1180   fvertexXY = new TH2F("fvertexXY","X-Y vertex position",30,0.,10.,30,0.,10.);
1181   fvertexZ = new TH1F("fvertexZ","Z vertex position",60,-30.,30.);
1182   fEvtMult = new TH1F("fEvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",100,0.,100.);
1183   fEvtMultvsJetPt = new TH2F("fEvtMultvsJetPt","Event multiplicity vs pT_{jet}",60,0.,60.,100,0.,100.);
1184   fPtvsEtaJet = new TH2F("fPtvsEtaJet","Pt vs #eta_{jet}",20,-1.,1.,60,0.,60.);
1185   fNpvsEtaJet = new TH2F("fNpvsEtaJet","N_{part} inside jet vs #eta_{jet}",20,-1.,1.,20,0,20);
1186   fNpevtvsEtaJet = new TH2F("fNpevtvsEtaJet","N_{part} in evt vs #eta_{jet}",20,-1.,1.,90,0,90);
1187   fPtvsPtJet = new TH2F("fPtvsPtJet","Pt vs #p_{Tjet}",60,0.,60.,60,0.,60.);
1188   fNpvsPtJet = new TH2F("fNpvsPtJet","N_{part} inside jet vs #pt_{Tjet}",60,0.,60.,20,0,20);
1189   fNpevtvsPtJet = new TH2F("fNpevtvsPtJet","N_{part} in evt vs #pt_{Tjet}",60,0.,60.,90,0,90);
1190   fPtvsPtJet1D = new TH1F("fPtvsPtJet1D","Pt vs #p_{Tjet}",60,0.,60.);
1191   fNpvsPtJet1D = new TH1F("fNpvsPtJet1D","N_{part} inside jet vs #pt_{Tjet}",60,0.,60.);
1192   fNpevtvsPtJet1D = new TH1F("fNpevtvsPtJet1D","N_{part} in evt vs #pt_{Tjet}",60,0.,60.);
1193   fptLeadingJet = new TH1F("fptLeadingJet","Pt leading Jet [GeV/c]",60,0.,60.);
1194   fetaLeadingJet = new TH1F("fetaLeadingJet","#eta leading jet",20,-1.,1.);
1195   fphiLeadingJet = new TH1F("fphiLeadingJet","#phi leading jet",12,0.,2*TMath::Pi());
1196   fptJet = new TH1F("fptJet","Pt Jets [GeV/c]",60,0.,60.);
1197   fetaJet = new TH1F("fetaJet","#eta jet",20,-1.,1.);
1198   fphiJet = new TH1F("fphiJet","#phi jet",12,0.,2*TMath::Pi());
1199   fNBadRunsH = new TH1F("fNBadRunsH","Number of events with Z vertex out of range", 1, 0., 1.);
1200
1201   SetProperties(fvertexXY, "vtx X", "vtx Y");
1202   SetProperties(fvertexZ, "vtx Z", "Count");
1203   SetProperties(fEvtMult, "N_{part} / event", "Count");
1204   SetProperties(fEvtMultvsJetPt, "p_{T jet}", "Event multiplicity");
1205   SetProperties(fPtvsEtaJet, "#eta_{leading jet}", "p_{T} part [GeV/c]");
1206   SetProperties(fNpvsEtaJet, "#eta_{leading jet}", "N_{part} in leading jet");
1207   SetProperties(fNpevtvsEtaJet, "#eta_{leading jet}", "N_{part} in event");
1208   SetProperties(fPtvsPtJet, "#p_{T leading jet}", "p_{T} part [GeV/c]");
1209   SetProperties(fNpvsPtJet, "#p_{T leading jet}", "N_{part} in leading jet");
1210   SetProperties(fNpevtvsPtJet, "#p_{T leading jet}", "N_{part} in event");
1211   SetProperties(fPtvsPtJet1D, "#p_{T leading jet}", "<p_{T}> part [GeV/c]");
1212   SetProperties(fNpvsPtJet1D, "#p_{T leading jet}", "<N_{part}> in leading jet");
1213   SetProperties(fNpevtvsPtJet1D, "#p_{T leading jet}", "<N_{part}> in event");
1214   SetProperties(fptLeadingJet, "p_{T} leading jet", "dN/dp_{T} leading jet");
1215   SetProperties(fetaLeadingJet, "#eta leading jet", "dN/d#eta leading jet");
1216   SetProperties(fphiLeadingJet, "#phi leading jet", "dN/d#phi leading jet");
1217   SetProperties(fptJet, "p_{T} jet [GeV/c]", "dN/dp_{T}");
1218   SetProperties(fetaJet, "#eta jet", "dN/d#eta");
1219   SetProperties(fphiJet, "#phi jet", "dN/d#phi");
1220
1221 }
1222
1223
1224 //////////////////////////////////////////////////////////////////////////////
1225 //////////////////////////////////////////////////////////////////////////////
1226
1227 void AliAnalysisTaskFragmentationFunction::FillMonoJetH(Int_t goodBin, AliAODJet* recJets, TClonesArray* tracks)
1228 {
1229   if (goodBin == 999) return;
1230
1231   Int_t nTracks = tracks->GetEntries();
1232   Float_t xi,t1,ene,dr2,deta2,dphi2, z, cosTheta, theta, kt;
1233   Bool_t jetOk1 = 0;
1234   Bool_t jetOk2 = 0;
1235   Bool_t jetOk3 = 0;
1236   Bool_t jetOk4 = 0;
1237   Bool_t jetOk5 = 0;
1238   Bool_t jetOk6 = 0;
1239   Bool_t jetOk7 = 0;
1240   Bool_t jetOk8 = 0;
1241   Bool_t jetOk9 = 0;
1242   Bool_t jetOk10 = 0;
1243
1244   fEtaMonoJet1H[goodBin]->Fill(recJets[0].Eta());
1245   fPhiMonoJet1H[goodBin]->Fill(recJets[0].Phi());
1246   fPtMonoJet1H[goodBin]->Fill(recJets[0].Pt());
1247   fEMonoJet1H[goodBin]->Fill(recJets[0].E());
1248
1249   Int_t mult = 0;
1250   for(Int_t it=0; it<nTracks; it++)
1251   {
1252     AliAODTrack* aodTrack = (AliAODTrack*)tracks->At(it);
1253     // Apply track cuts
1254     UInt_t status = aodTrack->GetStatus();
1255     if (status == 0) continue;
1256     if((fFilterMask>0)&&!(aodTrack->TestFilterBit(fFilterMask)))continue;
1257     mult++;
1258     Float_t etaT = aodTrack->Eta();
1259     Float_t phiT = aodTrack->Phi();
1260     Float_t ptT = aodTrack->Pt();
1261     // For Theta distribution
1262     Float_t pxT = aodTrack->Px();
1263     Float_t pyT = aodTrack->Py();
1264     Float_t pzT = aodTrack->Pz();
1265     Float_t pT = aodTrack->P();
1266     // Compute theta
1267     cosTheta = (pxT*recJets[0].Px()+pyT*recJets[0].Py()+pzT*recJets[0].Pz())/(pT*recJets[0].P());
1268     theta = TMath::ACos(cosTheta);
1269     // Compute xi
1270     deta2 = etaT - recJets[0].Eta();
1271     dphi2 = phiT - recJets[0].Phi();
1272     if (dphi2 < -TMath::Pi()) dphi2= -dphi2 - 2.0 * TMath::Pi();
1273     if (dphi2 > TMath::Pi()) dphi2 = 2.0 * TMath::Pi() - dphi2;
1274     dr2 = TMath::Sqrt(deta2 * deta2 + dphi2 * dphi2);
1275     t1  = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etaT)));
1276     ene = ptT*TMath::Sqrt(1.+1./(t1*t1));
1277     xi  = (Float_t) TMath::Log(recJets[0].E()/ptT);
1278     // Compute z
1279     z   = (Double_t)(ptT/recJets[0].E());
1280     // Compute kT/jT
1281     TVector3 partP; TVector3 jetP;
1282     jetP[0] = recJets[0].Px();
1283     jetP[1] = recJets[0].Py();
1284     jetP[2] = recJets[0].Pz();
1285     partP.SetPtEtaPhi(ptT,etaT,phiT);
1286     kt = TMath::Sin(partP.Angle(jetP))*partP.Mag();
1287     // Compute Jet shape
1288
1289     for(Int_t i2 = 0; i2 < fnRBin; i2++)
1290     {
1291       if ((dr2<farrayRadii[i2]) && (ptT > fPartPtCut))
1292       {
1293         if (i2 == 0) jetOk1 = 1;
1294         if (i2 == 1) jetOk2 = 1;
1295         if (i2 == 2) jetOk3 = 1;
1296         if (i2 == 3) jetOk4 = 1;
1297         if (i2 == 4) jetOk5 = 1;
1298         if (i2 == 5) jetOk6 = 1;
1299         if (i2 == 6) jetOk7 = 1;
1300         if (i2 == 7) jetOk8 = 1;
1301         if (i2 == 8) jetOk9 = 1;
1302         if (i2 == 9) jetOk10 = 1;
1303
1304         fdNdXiMonoJet1H[goodBin][i2]->Fill(xi);
1305         fdNdPtMonoJet1H[goodBin][i2]->Fill(ptT);
1306         fdNdZMonoJet1H[goodBin][i2]->Fill(z);
1307         fdNdThetaMonoJet1H[goodBin][i2]->Fill(theta);
1308         fdNdcosThetaMonoJet1H[goodBin][i2]->Fill(cosTheta);
1309         fdNdkTMonoJet1H[goodBin][i2]->Fill(kt);
1310         fdNdpTvsZMonoJet1H[goodBin][i2]->Fill(z,1/((fPthadHBinMax-fPthadHBinMin)/fnPthadHBin));
1311
1312         fThetaPtPartMonoJet1H[goodBin][i2]->Fill(ptT,theta);
1313         fcosThetaPtPartMonoJet1H[goodBin][i2]->Fill(ptT,cosTheta);
1314         fkTPtPartMonoJet1H[goodBin][i2]->Fill(ptT,kt);
1315         fThetaPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),theta);
1316         fcosThetaPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),cosTheta);
1317         fkTPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),kt);
1318         fpTPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),ptT);
1319       }
1320     }
1321   }
1322   fEvtMultvsJetPt->Fill(recJets[0].Pt(),mult);
1323
1324   if (jetOk1)  fNMonoJet1sH[goodBin][0]->Fill(0.5);
1325   if (jetOk2)  fNMonoJet1sH[goodBin][1]->Fill(0.5);
1326   if (jetOk3)  fNMonoJet1sH[goodBin][2]->Fill(0.5);
1327   if (jetOk4)  fNMonoJet1sH[goodBin][3]->Fill(0.5);
1328   if (jetOk5)  fNMonoJet1sH[goodBin][4]->Fill(0.5);
1329   if (jetOk6)  fNMonoJet1sH[goodBin][5]->Fill(0.5);
1330   if (jetOk7)  fNMonoJet1sH[goodBin][6]->Fill(0.5);
1331   if (jetOk8)  fNMonoJet1sH[goodBin][7]->Fill(0.5);
1332   if (jetOk9)  fNMonoJet1sH[goodBin][8]->Fill(0.5);
1333   if (jetOk10) fNMonoJet1sH[goodBin][9]->Fill(0.5);
1334
1335 }
1336
1337 //////////////////////////////////////////////////////////////////////////////
1338 //////////////////////////////////////////////////////////////////////////////
1339
1340 void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y)
1341 {
1342   //Set properties of histos (x and y title and error propagation)
1343   h->SetXTitle(x);
1344   h->SetYTitle(y);
1345   h->GetXaxis()->SetTitleColor(1);
1346   h->GetYaxis()->SetTitleColor(1);
1347   h->Sumw2();
1348 }
1349
1350 //////////////////////////////////////////////////////////////////////////////
1351 //////////////////////////////////////////////////////////////////////////////
1352
1353 void AliAnalysisTaskFragmentationFunction::MakeJetContainer()
1354 {
1355   //
1356   // Create the particle container for the correction framework manager and
1357   // link it
1358   //
1359   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1360   const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
1361   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1362   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1363
1364   // can we neglect migration in eta and phi?
1365   // phi should be no problem since we cover full phi and are phi symmetric
1366   // eta migration is more difficult  due to needed acceptance correction
1367   // in limited eta range
1368
1369   //arrays for the number of bins in each dimension
1370   Int_t iBin[kNvar];
1371   iBin[0] = 100; //bins in pt
1372   iBin[1] =  1; //bins in eta
1373   iBin[2] = 1; // bins in phi
1374
1375   //arrays for lower bounds :
1376   Double_t* binEdges[kNvar];
1377   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1378     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1379
1380   //values for bin lower bounds
1381   //  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);
1382   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1383   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1384   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1385
1386   /*
1387   for(int i = 0;i < kMaxStep*2;++i){
1388     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
1389     for (int k=0; k<kNvar; k++) {
1390       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1391     }
1392   }
1393   //create correlation matrix for unfolding
1394   Int_t thnDim[2*kNvar];
1395   for (int k=0; k<kNvar; k++) {
1396     //first half  : reconstructed
1397     //second half : MC
1398     thnDim[k]      = iBin[k];
1399     thnDim[k+kNvar] = iBin[k];
1400   }
1401
1402   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1403   for (int k=0; k<kNvar; k++) {
1404     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1405     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1406   }
1407   fhnCorrelation->Sumw2();
1408
1409   // Add a histogram for Fake jets
1410   //  thnDim[3] = AliPID::kSPECIES;
1411   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
1412   // for(Int_t idim = 0; idim < kNvar; idim++)
1413   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
1414   */
1415 }
1416
1417 //////////////////////////////////////////////////////////////////////////////
1418 //////////////////////////////////////////////////////////////////////////////