PWGJE/UserTasks/AliAnalysisTaskJetShape.cxx

   1 // Class AliAnalysisTaskJetShape
   2 // Author martin.poghosyan@cern.ch
   3 //
   4
   5 #include "TChain.h"
   6 #include "TTree.h"
   7 #include "TH1.h"
   8 #include "TH1F.h"
   9 #include "TH2F.h"
  10 #include "TH3F.h"
  11 #include "TCanvas.h"
  12 #include "TParticlePDG.h"
  13 #include "TParticle.h"
  14 #include "TRandom.h"
  15 #include "TDatabasePDG.h"
  16 #include <TPDGCode.h>
  17 #include "THnSparse.h"
  18 #include "TGraph.h"
  19 #include "TArrayF.h"
  20 #include "TArrayD.h"
  21 #include <TVector3.h>
  22 #include <TClonesArray.h>
  23 #include "AliLog.h"
  24
  25 #include "AliCDBManager.h"
  26 #include "AliGeomManager.h"
  27 #include "AliCDBEntry.h"
  28 #include "AliCDBPath.h"
  29 #include "AliAlignObjParams.h"
  30 #include "AliAnalysisTaskSE.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliESDEvent.h"
  33 #include "AliESDInputHandler.h"
  34 #include "AliMultiplicity.h"
  35 #include "AliESDFMD.h"
  36 #include "AliESDVZERO.h"
  37 #include "AliESDInputHandler.h"
  38 #include "AliVEvent.h"
  39 #include "AliCentrality.h"
  40 #include "AliAnalysisHelperJetTasks.h"
  41 #include "AliInputEventHandler.h"
  42 #include "AliAODJetEventBackground.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODHandler.h"
  45 #include "AliAODJet.h"
  46 #include "AliAODEvent.h"
  47 #include "AliAODTrack.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliGenEventHeader.h"
  50 #include "AliHeader.h"
  51 #include "AliStack.h"
  52 #include "AliMCEvent.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliInputEventHandler.h"
  55 #include "AliLHCData.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliGenDPMjetEventHeader.h"
  58 #include "AliESDtrackCuts.h"
  59 #include "AliTriggerAnalysis.h"
  60
  61
  62
  63 #include "AliAnalysisTaskJetShape.h"
  64
  65 /////////////////////////////////////
  66
  67 ClassImp(AliAnalysisTaskJetShape)
  68
  69
  70   AliAnalysisTaskJetShape::AliAnalysisTaskJetShape(const char *name)
  71 : AliAnalysisTaskSE(name),
  72   //  fOutputTree(0),
  73   fOutputList(0),
  74   fESD(0x0),
  75   fAODIn(0x0),
  76   fAODOut(0x0),
  77   fAODExtension(0x0),
  78   fkMC(0),
  79   fCMSE(0),
  80   fRunNb(0),
  81   fkIsPhysSel(0),
  82   fNonStdFile(""),
  83   fFilterMask(0),
  84   fOfflineTrgMask(AliVEvent::kAny),
  85   fCentMin(0.),
  86   fCentMax(80.),
  87   fEvtClassMin(0),
  88   fEvtClassMax(4),
  89   fPtJcorrMin(20),
  90   fPtJBcorrMin(20),
  91   fJpPtmin(1),
  92   fJaPtmin(0.5),
  93   fVtxMinContrib(1),
  94   fVtxZMin(-10),
  95   fVtxZMax(10),
  96   fkIsPbPb(0),
  97   fhPtJL(0),
  98   fhAreaJL(0),
  99   fanalJetSubStrHM(0),
 100   fanalJetSubStrMCHM(0)
 101 {
 102
 103
 104   fgkbinCent[0] = -0.099;
 105   fgkbinCent[1] = 5;
 106   fgkbinCent[2] = 10;
 107   fgkbinCent[3] = 20;
 108   fgkbinCent[4] = 40;
 109   fgkbinCent[5] = 60;
 110   fgkbinCent[6] = 80;
 111   fgkbinCent[7] =100;
 112
 113   fBackgroundBranch[0] = "",
 114   fBackgroundBranch[1] = "",
 115   fJetBranchName[0] = "";
 116   fJetBranchName[1] = "";
 117
 118   for(Int_t i=0; i<3; i++)
 119     {
 120       fhPtresVsPt[i]=0;
 121       fhPhiresVsPhi[i]=0;
 122       fhEtaresVsEta[i]=0;
 123       fhDCAxy[i]=0;
 124       fhDCAz[i]=0;
 125       fhTrackPtEtaPhi[i]=0;
 126     }
 127   //  fkIsPbPb = kFALSE;
 128   //  fDebug=0;
 129
 130
 131
 132   fanalJetSubStrHM   = new AliAnalysisTaskJetShapeHM("rec", kFALSE);
 133   fanalJetSubStrMCHM = new AliAnalysisTaskJetShapeHM("truth", kTRUE);
 134
 135   DefineOutput(1, TList::Class());
 136
 137   //   DefineOutput(1, TTree::Class());
 138 }
 139
 140
 141 AliAnalysisTaskJetShape::~AliAnalysisTaskJetShape()
 142 {
 143    if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 144     printf("Deleteing output\n");
 145
 146     if(fOutputList){
 147       delete fOutputList;
 148       fOutputList = 0;
 149     }
 150
 151     //    if(fTriggerAnalysis)
 152     //      delete fTriggerAnalysis;
 153
 154    }
 155 }
 156
 157
 158 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool():TObject(),
 159   fvecJ(0,0,0),
 160   fvecB1(0,0,0),
 161   fvecB2(0,0,0),
 162   fRmax(0),
 163   fPRecInRJ(0,0,0),
 164   fList(0),
 165   fEntries(0)
 166 {
 167   fList = new TClonesArray("TVector3",10000);
 168
 169   TVector3 v(0,0,0);
 170   SetVecJ(v);
 171   fRmax = -0.5;
 172   fR[0] =0.1;
 173   fR[1] =0.2;
 174   fR[2] =0.3;
 175   fEntries=0;
 176
 177   for(Int_t i1=0; i1<fgkbinR; i1++)
 178     {
 179       for(Int_t i2=0; i2<2; i2++)
 180         {
 181           fListJ[i1][i2].Set(1000);
 182           fListB1[i1][i2].Set(1000);
 183           fListB2[i1][i2].Set(1000);
 184           fListJc[i1][i2]  = 0;
 185           fListB1c[i1][i2] = 0;
 186           fListB2c[i1][i2] = 0;
 187           //      fListJ[i1][i2]("TVector3",10000);
 188           //      fListB1[i1][i2] = TClonesArray("TVector3",10000);
 189           //      fListB2[i1][i2] = TClonesArray("TVector3",10000);
 190
 191                   fListJ[i1][i2].Reset(-1);
 192                   fListB1[i1][i2].Reset(-1);
 193                   fListB2[i1][i2].Reset(-1);
 194
 195                   fPRecJ[i1][i2].SetXYZ(0,0,0);
 196
 197         }
 198     }
 199
 200   fPRecInRJ.SetXYZ(0,0,0);
 201
 202   fPtMinTr[0] = 2;
 203   fPtMinTr[1] = 0.5;
 204
 205
 206 }
 207
 208 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool(TClonesArray *list):TObject(),
 209   fvecJ(0,0,0),
 210   fvecB1(0,0,0),
 211   fvecB2(0,0,0),
 212   fRmax(0),
 213   fPRecInRJ(0,0,0),
 214   fList(list),
 215   fEntries(0)
 216 {
 217
 218   TVector3 v(0,0,0);
 219   SetVecJ(v);
 220   fRmax = 0.4;
 221   fR[0] =0.1;
 222   fR[1] =0.2;
 223   fR[2] =0.3;
 224   fEntries=0;
 225
 226   for(Int_t i1=0; i1<fgkbinR; i1++)
 227     {
 228       for(Int_t i2=0; i2<2; i2++)
 229         {
 230           fListJ[i1][i2].Set(1000);
 231           fListB1[i1][i2].Set(1000);
 232           fListB2[i1][i2].Set(1000);
 233           fListJc[i1][i2]  = 0;
 234           fListB1c[i1][i2] = 0;
 235           fListB2c[i1][i2] = 0;
 236           //      fListJ[i1][i2]("TVector3",10000);
 237           //      fListB1[i1][i2] = TClonesArray("TVector3",10000);
 238           //      fListB2[i1][i2] = TClonesArray("TVector3",10000);
 239                   fListJ[i1][i2].Reset(-1);
 240                   fListB1[i1][i2].Reset(-1);
 241                   fListB2[i1][i2].Reset(-1);
 242                   fPRecJ[i1][i2].SetXYZ(0,0,0);
 243
 244         }
 245     }
 246
 247   fPRecInRJ.SetXYZ(0,0,0);
 248
 249   fPtMinTr[0] = 2;
 250   fPtMinTr[1] = 0.5;
 251
 252 }
 253
 254
 255
 256
 257 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::~AliAnalysisTaskJetShapeTool()
 258 {
 259   fList->Delete();
 260   fEntries=0;
 261
 262   //  if(fList)
 263   //    delete fList;
 264
 265 }
 266
 267
 268 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::SetVecJ(TVector3 v)
 269 {
 270   fvecJ.SetXYZ(v.X(), v.Y(), v.Z());
 271   fvecB1.SetXYZ(-v.Y(), v.X(), v.Z());
 272   fvecB2.SetXYZ(v.Y(),-v.X(), v.Z());
 273 }
 274
 275
 276 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Clean()
 277 {
 278   //  printf("AnalJetSubStrTool::Clean()\n");
 279   /*
 280   fList->Delete();
 281   for(Int_t i1=0; i1<fgkbinR; i1++)
 282     {
 283       for(Int_t i2=0; i2<2; i2++)
 284         {
 285           fListJ[i1][i2]->Delete();
 286           fListB1[i1][i2]->Delete();
 287           fListB2[i1][i2]->Delete();
 288         }
 289     }
 290
 291   */
 292
 293   //  fList.SetOwner();
 294   //  fList->Delete();
 295   //  fEntries=0;
 296   fPRecInRJ.SetXYZ(0,0,0);
 297
 298   for(Int_t i1=0; i1<fgkbinR; i1++)
 299     {
 300       for(Int_t i2=0; i2<2; i2++)
 301         {
 302                   fListJ[i1][i2].Reset(-1);
 303                   fListB1[i1][i2].Reset(-1);
 304                   fListB2[i1][i2].Reset(-1);
 305           //      fListJ[i1][i2].Set(0);
 306           //      fListB1[i1][i2].Set(0);
 307           //      fListB2[i1][i2].Set(0);
 308           fListJc[i1][i2]  = 0;
 309           fListB1c[i1][i2] = 0;
 310           fListB2c[i1][i2] = 0;
 311
 312           fPRecJ[i1][i2].SetXYZ(0,0,0);
 313         }
 314     }
 315
 316 }
 317
 318
 319 //void AnalJetSubStrTool::Print(Option_t* /*option*/) const
 320 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT() const
 321 {
 322
 323   //      PRINT(fList, "all");
 324
 325           //  fList->Print();
 326   for(Int_t i1=0; i1<fgkbinR; i1++)
 327     {
 328
 329       if(i1!=0) continue;
 330      for(Int_t i2=0; i2<2; i2++)
 331         {
 332
 333           //      printf("# %d   %d\n",i1, i2);
 334           PRINT(fListJ[i1][i2], fListJc[i1][i2], Form("J_%d_%d",i1,i2));
 335                   //              PRINT(fListB1[i1][i2], Form("B1_%d_%d",i1,i2));
 336                   //              PRINT(fListB2[i1][i2], Form("B2_%d_%d",i1,i2));
 337           //      fListJ[i1][i2].Print("",1);
 338           //      fListB1[i1][i2]->Print();
 339           //      fListB2[i1][i2]->Print();
 340         }
 341     }
 342
 343 }
 344
 345 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT(TArrayI a, Int_t n, const char *txt) const
 346 {
 347   printf("%s :%d \n",txt, n);
 348   Int_t count = 0;
 349   for(Int_t i=0; i<n; i++){
 350     printf("%4d   ", a.At(i));
 351
 352     if(count==20)
 353       {
 354         printf("\n");
 355         count=0;
 356       }
 357     else count++;
 358   }
 359   if(count!=20) printf("\n");
 360
 361 }
 362
 363
 364
 365  void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Init()
 366 {
 367   Int_t Nev = fEntries;
 368
 369   //    PRINT();
 370
 371   //  printf("Nev = %d\n",Nev);
 372
 373   for(Int_t iP=0; iP<Nev; iP++)
 374     {
 375       TVector3 *v = (TVector3*) fList->At(iP);
 376       if(!v) {
 377         printf("ERROR : entry #%d not found\n",iP);
 378         continue;
 379       }
 380
 381       //      printf("#%3d   ",iP); v->Print();
 382
 383       Double_t R = CalcR(fvecJ, *v);
 384       //     printf("R = %f\n",R);
 385       if(R<fRmax)
 386         {
 387           for(Int_t iT = 0; iT < fgkbinR; iT++)
 388             {
 389               Int_t type = 0;
 390               if(R>fR[iT]) type = 1;
 391
 392               if(v->Pt() < fPtMinTr[type]) continue;
 393               fPRecJ[iT][type]+=*v;
 394               AddToJ(iT, type, iP);
 395             }
 396
 397           if(v->Pt() < fPtMinTr[0]) continue;
 398           fPRecInRJ+=*v;
 399
 400           continue;
 401         }
 402
 403       R = CalcR(fvecB1, *v);
 404       if(R<fRmax)
 405         {
 406           for(Int_t iT = 0; iT < fgkbinR; iT++)
 407             {
 408               Int_t type = 0;
 409               if(R>fR[iT]) type = 1;
 410
 411               if(v->Pt() < fPtMinTr[type]) continue;
 412
 413               AddToB1(iT, type, iP);
 414             }
 415           continue;
 416         }
 417
 418       R = CalcR(fvecB2, *v);
 419       if(R<fRmax)
 420         {
 421           for(Int_t iT = 0; iT < fgkbinR; iT++)
 422             {
 423               Int_t type = 0;
 424               if(R>fR[iT]) type = 1;
 425
 426               if(v->Pt() < fPtMinTr[type]) continue;
 427
 428               AddToB2(iT, type, iP);
 429             }
 430           continue;
 431         }
 432     }
 433
 434
 435
 436   /*
 437   for(Int_t i1=0; i1<fgkbinR; i1++)
 438     {
 439       for(Int_t i2=0; i2<2; i2++)
 440         {
 441           //      if(fListJc[i1][i2]<2) fListJc[i1][i2]=0;
 442           //      if(fListB1c[i1][i2]<2) fListB1c[i1][i2]=0;
 443           //      if(fListB2c[i1][i2]<2) fListB2c[i1][i2]=0;
 444           fListJ[i1][i2].Set(fListJc[i1][i2]);
 445           fListB1[i1][i2].Set(fListB1c[i1][i2]);
 446           fListB2[i1][i2].Set(fListB2c[i1][i2]);
 447         }
 448     }
 449   */
 450
 451 }
 452
 453
 454
 455 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcR(TVector3 v1, TVector3 v2)
 456 {
 457
 458   Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
 459   //  dphi*=(0.9/TMath::Pi());
 460   Double_t deta = v1.Eta() - v2.Eta();
 461   Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
 462   return RB;
 463 }
 464
 465
 466 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcdPhi(Double_t phi1, Double_t phi2)
 467 {
 468
 469   while(phi1<0) phi1+=TMath::TwoPi();
 470   while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
 471
 472   while(phi2<0) phi2+=TMath::TwoPi();
 473   while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
 474
 475   Double_t dphi = phi1- phi2;
 476   if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
 477   if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
 478
 479   return dphi;
 480 }
 481
 482
 483 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::GetLocalMom(TVector3 vecJ, TVector3 *q)
 484 {
 485   // theta and phi of a particle in loc.syst of fvecJ
 486
 487
 488   Double_t PT = vecJ.Perp();
 489   Double_t P0 = vecJ.Mag();
 490
 491   Double_t q0=q->X();
 492   Double_t q1=q->Y();
 493   Double_t q2=q->Z();
 494
 495   Double_t qx1 = vecJ(2)/P0/PT*(vecJ(0)*q0+vecJ(1)*q1) - PT/P0*q2;
 496   Double_t qy1 = -vecJ(1)/PT*q0 + vecJ(0)/PT*q1;
 497   Double_t qz1 = 1/P0*(vecJ(0)*q0+vecJ(1)*q1+vecJ(2)*q2);
 498
 499   //  Double_t q00 = TMath::Sqrt(qx1*qx1 + qy1*qy1 + qz1*qz1);
 500   //  phi = TMath::ATan2(qy1, qx1);
 501   //  if(phi<0) phi+=fTwoPi;
 502   //  theta = TMath::ACos(qz1/q00);
 503
 504   q->SetXYZ(qx1, qy1, qz1);
 505
 506   return;
 507
 508 }
 509
 510 /*
 511
 512 Bool_t AnalJetSubStrTool::FindCorrelationAxesAnd(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
 513 {
 514   Double_t TwoPi = TMath::TwoPi();
 515
 516 //
 517 // 1st track loop to determine the sphericity matrix
 518 //
 519       // Initialize Shericity matrix
 520       Float_t mxx    = 0.;
 521       Float_t myy    = 0.;
 522       Float_t mxy    = 0.;
 523       Int_t   nc     = 0;
 524       Float_t sump2  = 0.;
 525       Float_t ptMax  = 0.;
 526       Float_t etaMax = 0.;
 527       Float_t phiMax = 0.;
 528       Int_t   iMax   = -1;
 529       Float_t ptsum  = 0.;
 530
 531
 532   Int_t Nev =  list.GetSize();
 533   if( Nev <2) return kFALSE;
 534
 535   Int_t indexpmax = -1;
 536   Double_t pmax = -1;
 537   Double_t phipmax = 0;
 538
 539   Int_t indexpzmax = -1;
 540   Double_t pzmax = -1;
 541   Double_t phipzmax = 0;
 542
 543
 544   Int_t indexthetamin = -1;
 545   Double_t thetamin = 1000;
 546   Double_t phithetamin = 0;
 547
 548   for(Int_t ip=0; ip< Nev; ip++) {
 549
 550     TVector3* part = (TVector3*) GetAt(list.At(ip));
 551     if(!part) continue;
 552
 553     TVector3 vtmp(*part);
 554     GetLocalMom(vec, &vtmp);
 555
 556     Float_t ppjX = vtmp.X();
 557     Float_t ppjY = vtmp.Y();
 558     Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
 559
 560
 561     Float_t pt = ppjT;//part->Z();
 562     Float_t deta = part->Eta();
 563     Float_t dphi = part->Phi();
 564
 565               mxx += (ppjX * ppjX / ppjT);
 566               myy += (ppjY * ppjY / ppjT);
 567               mxy += (ppjX * ppjY / ppjT);
 568               nc++;
 569               sump2 += ppjT;
 570
 571
 572     if(vtmp.Mag()>pmax)
 573       {
 574         pmax=vtmp.Mag();
 575         indexpmax = ip;
 576         phipmax=vtmp.Phi();
 577       }
 578
 579     if(vtmp.Z()>pzmax)
 580       {
 581         pzmax=vtmp.Z();
 582         indexpzmax = ip;
 583         phipzmax=vtmp.Phi();
 584       }
 585     if(vtmp.Theta()<thetamin)
 586       {
 587         thetamin=vtmp.Theta();
 588         indexthetamin = ip;
 589         phithetamin=vtmp.Phi();
 590       }
 591
 592
 593           }
 594
 595
 596
 597
 598 //
 599 // At this point we have mxx, myy, mxy
 600   if (nc == 0) return kFALSE;
 601 // Shericity Matrix
 602       const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
 603       TMatrixDSym m0(2, ele);
 604 // Find eigenvectors
 605       TMatrixDSymEigen m(m0);
 606       TVectorD eval(2);
 607       TMatrixD evecm = m.GetEigenVectors();
 608       eval  = m.GetEigenValues();
 609 // Largest eigenvector
 610       Int_t jev = 0;
 611       if (eval[0] < eval[1]) jev = 1;
 612       TVectorD evec0(2);
 613 // Principle axis
 614       evec0 = TMatrixDColumn(evecm, jev);
 615       TVector2 evec(evec0[0], evec0[1]);
 616 // Principle axis from leading partice
 617
 618 //      Float_t phiM = TMath::ATan2(phiMax, etaMax);
 619 //      TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
 620 //      Float_t phistM = evecM.DeltaPhi(evec);
 621 //      if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
 622
 623
 624       Double_t phiTA = evec.Phi();
 625
 626
 627       Double_t phiDir = 0;
 628   if(scenario ==0)
 629     {
 630       phiDir = phipmax;
 631     }
 632   else if(scenario ==1)
 633     {
 634       phiDir = phipzmax;
 635     }
 636
 637   else if(scenario ==2)
 638     {
 639       phiDir = phithetamin;
 640     }
 641   else
 642     return kFALSE;
 643
 644
 645       Phi = phiTA;
 646
 647       if( TMath::Abs(CalcdPhi(phiDir, phiTA)) <TMath::Pi()/2)
 648         return kTRUE;
 649       else {
 650         Phi+=TMath::Pi();
 651         if(Phi>TwoPi) Phi-=TwoPi;
 652         return kTRUE;
 653       }
 654
 655
 656  return kTRUE;
 657 }
 658 */
 659
 660 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxes(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
 661 {
 662   // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
 663
 664   Double_t xmean , ymean;
 665
 666   Double_t TwoPi = TMath::TwoPi();
 667
 668   xmean = 0;
 669   ymean = 0;
 670   Phi= -999;
 671   //find axes
 672
 673   Double_t x=0;
 674   Double_t x2=0;
 675   Double_t y=0;
 676   Double_t y2=0;
 677   Double_t xy=0;
 678
 679   //  Double_t phiLoc, thetaLoc;
 680
 681   Int_t N=0;
 682
 683   Int_t Nev =  list.GetSize();
 684   if( Nev <2) return kFALSE;
 685
 686
 687
 688
 689   Int_t Index = -1;
 690   Double_t val = -1;
 691   if(scenario==2) val = 100;
 692   Double_t phiDir = -99;
 693
 694
 695
 696   for(Int_t ip=0; ip< Nev; ip++) {
 697     if(list.At(ip)<0) break;
 698
 699     TVector3* part = (TVector3*) GetAt(list.At(ip));
 700     if(!part) continue;
 701
 702
 703     Double_t xx, yy;
 704     /*
 705     TVector3 vtmp(*part);
 706     GetLocalMom(vec, &vtmp);
 707     Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
 708     Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
 709     */
 710
 711     xx = CalcdPhi(part->Phi(), vec.Phi());
 712     yy = part->Eta() - vec.Eta();
 713
 714     x+=xx;
 715     y+=yy;
 716     x2+=(xx*xx);
 717     y2+=(yy*yy);
 718     xy+=(xx*yy);
 719     N++;
 720
 721     /*
 722   if(scenario ==0)
 723     {
 724       if(vtmp.Mag()>val)
 725         {
 726           val=vtmp.Mag();
 727           Index = ip;
 728           phiDir=vtmp.Phi();
 729         }
 730     }
 731   else if(scenario ==1)
 732     {
 733       if(vtmp.Z()>val)
 734         {
 735           val=vtmp.Z();
 736           Index = ip;
 737           phiDir=vtmp.Phi();
 738         }
 739     }
 740   else if(scenario ==2)
 741     {
 742       if(vtmp.Theta()<val)
 743         {
 744           val=vtmp.Theta();
 745           Index = ip;
 746           phiDir=vtmp.Phi();
 747         }
 748     }
 749   else
 750 */
 751 if(scenario ==3)
 752     {
 753       if(part->Pt()> val)
 754         {
 755           val =part->Pt();
 756           Index = ip;
 757           //      phiDir=vtmp.Phi();
 758           phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
 759         }
 760     }
 761
 762
 763   }
 764
 765   if(N<2) return kFALSE;
 766
 767   //  return kFALSE;
 768
 769   x/=N;
 770   y/=N;
 771   x2/=N;
 772   y2/=N;
 773   xy/=N;
 774
 775
 776   Double_t A = xy-x*y;
 777   Double_t B = x2-x*x+y*y-y2;
 778   Double_t D = TMath::Sqrt(B*B+4*A*A);
 779
 780   //  printf("N = %d  A = %f\n",N, A);
 781   //  list.Print();
 782
 783   Double_t a1 = (-B+D)/2/A;
 784   //  Double_t a2 = (-B-D)/2/A;
 785   //  Double_t b1 = y-a1*x;
 786   //  Double_t b2 = y-a2*x;
 787   //  Double_t phiDir = TMath::ATan2(y, x);
 788   // while(phiDir<0) phiDir+=TwoPi;
 789     Double_t phi = TMath::ATan(a1);
 790
 791   while(phi>TwoPi) phi-=TwoPi;
 792   while(phi< 0 )   phi+=TwoPi;
 793
 794   Phi=CalcdPhi0To2pi(phi, 0);
 795
 796   if(scenario!=4 && Index<0) return kFALSE;
 797
 798   if(Index>=0)
 799     {
 800       if( TMath::Abs(CalcdPhi(phiDir, phi)) <TMath::Pi()/2)
 801         return kTRUE;
 802       else
 803         {
 804           phi+=TMath::Pi();
 805           Phi = CalcdPhi0To2pi(phi);
 806           return kTRUE;
 807         }
 808     }
 809
 810
 811   Double_t xmax = -100;
 812   Double_t xmin =  100;
 813
 814   xmean = x;
 815   ymean = y;
 816
 817
 818   for(Int_t ip=0; ip< Nev; ip++) {
 819     if(list.At(ip)<0) break;
 820
 821     TVector3* part = (TVector3*) GetAt(list.At(ip));
 822     if(!part) continue;
 823
 824     TVector3 vtmp(*part);
 825     GetLocalMom(vec, &vtmp);
 826
 827     Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
 828     Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
 829
 830     Double_t x1 = (xx-xmean)*TMath::Cos(Phi) +(yy-ymean)*TMath::Sin(Phi);
 831     if(x1 > xmax) xmax = x1;
 832     if(x1 < xmin) xmin = x1;
 833     //    printf("c3\n");
 834   }
 835
 836   if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
 837
 838
 839
 840
 841    while(Phi>TwoPi) Phi-=TwoPi;
 842
 843   return kTRUE;
 844 }
 845
 846
 847
 848
 849 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxesCorr(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario, Double_t R)
 850 {
 851   // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
 852
 853   Double_t xmean , ymean;
 854
 855   Double_t TwoPi = TMath::TwoPi();
 856
 857   xmean = 0;
 858   ymean = 0;
 859   Phi= -999;
 860   //find axes
 861
 862   Double_t x=0;
 863   Double_t x2=0;
 864   Double_t y=0;
 865   Double_t y2=0;
 866   Double_t xy=0;
 867
 868   //  Double_t phiLoc, thetaLoc;
 869
 870   Int_t N=0;
 871
 872   Int_t Nev =  list.GetSize();
 873   if( Nev <2) return kFALSE;
 874
 875
 876
 877
 878   Int_t Index = -1;
 879   Double_t val = -1;
 880   if(scenario==2) val = 100;
 881   Double_t phiDir = -99;
 882
 883
 884
 885   for(Int_t ip=0; ip< Nev; ip++) {
 886     if(list.At(ip)<0) break;
 887
 888     TVector3* part = (TVector3*) GetAt(list.At(ip));
 889     if(!part) continue;
 890
 891     Double_t xx, yy;
 892     /*
 893     TVector3 vtmp(*part);
 894     GetLocalMom(vec, &vtmp);
 895     Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
 896     Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
 897     */
 898
 899     xx = CalcdPhi(part->Phi(), vec.Phi());
 900     yy = part->Eta() - vec.Eta();
 901
 902     x+=xx;
 903     y+=yy;
 904     x2+=(xx*xx);
 905     y2+=(yy*yy);
 906     xy+=(xx*yy);
 907     N++;
 908
 909     /*
 910   if(scenario ==0)
 911     {
 912       if(vtmp.Mag()>val)
 913         {
 914           val=vtmp.Mag();
 915           Index = ip;
 916           phiDir=vtmp.Phi();
 917         }
 918     }
 919   else if(scenario ==1)
 920     {
 921       if(vtmp.Z()>val)
 922         {
 923           val=vtmp.Z();
 924           Index = ip;
 925           phiDir=vtmp.Phi();
 926         }
 927     }
 928   else if(scenario ==2)
 929     {
 930       if(vtmp.Theta()<val)
 931         {
 932           val=vtmp.Theta();
 933           Index = ip;
 934           phiDir=vtmp.Phi();
 935         }
 936     }
 937   else
 938 */
 939 if(scenario ==3)
 940     {
 941       if(part->Pt()> val)
 942         {
 943           val =part->Pt();
 944           Index = ip;
 945           //      phiDir=vtmp.Phi();
 946           phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
 947         }
 948     }
 949
 950
 951   }
 952
 953   if(N<2) return kFALSE;
 954   if(scenario!=4 && Index<0) return kFALSE;
 955
 956   //  return kFALSE;
 957
 958   x/=N;
 959   y/=N;
 960   x2/=N;
 961   y2/=N;
 962   xy/=N;
 963
 964
 965   Double_t A = xy-x*y;
 966   Double_t B = x2-x*x+y*y-y2;
 967   Double_t D = TMath::Sqrt(B*B+4*A*A);
 968
 969   //  printf("N = %d  A = %f\n",N, A);
 970   //  list.Print();
 971
 972   Double_t a1 = (-B+D)/2/A;
 973   //  Double_t a2 = (-B-D)/2/A;
 974   //  Double_t b1 = y-a1*x;
 975   //  Double_t b2 = y-a2*x;
 976   //  Double_t phiDir = TMath::ATan2(y, x);
 977   // while(phiDir<0) phiDir+=TwoPi;
 978     Double_t phi = TMath::ATan(a1);
 979
 980   Phi=CalcdPhi0To2pi(phi, 0);
 981   if( TMath::Abs(CalcdPhi(phiDir, phi))  > TMath::Pi()/2)
 982     Phi=CalcdPhi0To2pi(phi+TMath::Pi(), 0);
 983
 984
 985
 986
 987   if(Index>=0)
 988     {
 989       if(N==2)
 990         return kTRUE;
 991
 992
 993
 994       Double_t phiMinus = 0;
 995       Int_t Nminus = 0;
 996       for(Int_t ip=0; ip< Nev; ip++)
 997         {
 998           if(list.At(ip)<0) break;
 999
1000           TVector3* part = (TVector3*) GetAt(list.At(ip));
1001           if(!part) continue;
1002           Double_t xx0 = CalcdPhi(part->Phi(), vec.Phi());
1003           Double_t yy0 = part->Eta() - vec.Eta();
1004
1005           Double_t xx  = (x*N - xx0)/(N-1);
1006           Double_t yy  = (y*N - yy0)/(N-1);
1007           Double_t xx2 = (x2*N - xx0*xx0)/(N-1);
1008           Double_t yy2 = (y2*N - yy0*yy0)/(N-1);
1009           Double_t xxyy= (xy*N - xx0*yy0)/(N-1);
1010
1011
1012           Double_t AA = xxyy-xx*yy;
1013           Double_t BB = xx2-xx*xx+yy*yy-yy2;
1014           Double_t DD = TMath::Sqrt(BB*BB+4*AA*AA);
1015           Double_t aa1 = (-BB+DD)/2/AA;
1016           Double_t phi1 = TMath::ATan(aa1);
1017
1018           if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
1019             phi1+=TMath::Pi();
1020
1021           phiMinus+=CalcdPhi(Phi, phi1);
1022           Nminus++;
1023
1024         }
1025
1026
1027       Double_t phiPlus = 0;
1028       Int_t Nplus = 0;
1029       for(Int_t ip=0; ip< -Nev; ip++)
1030         {
1031           if(list.At(ip)<0) break;
1032
1033           TVector3* part = (TVector3*) GetAt(list.At(ip));
1034           if(!part) continue;
1035           Double_t rtmp = gRandom->Uniform(0, R);
1036           Double_t phitmp = gRandom->Uniform(0, TMath::TwoPi());
1037           Double_t xx0 = rtmp*TMath::Cos(phitmp);
1038           Double_t yy0 = rtmp*TMath::Sin(phitmp);
1039
1040           Double_t xx  = (x*N + xx0)/(N+1);
1041           Double_t yy  = (y*N + yy0)/(N+1);
1042           Double_t xx2 = (x2*N + xx0*xx0)/(N+1);
1043           Double_t yy2 = (y2*N + yy0*yy0)/(N+1);
1044           Double_t xxyy= (xy*N + xx0*yy0)/(N+1);
1045
1046
1047           Double_t AA = xxyy-xx*yy;
1048           Double_t BB = xx2-xx*xx+yy*yy-yy2;
1049           Double_t DD = TMath::Sqrt(BB*BB+4*AA*AA);
1050           Double_t aa1 = (-BB+DD)/2/AA;
1051           Double_t phi1 = TMath::ATan(aa1);
1052
1053           if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
1054             phi1+=TMath::Pi();
1055
1056           phiPlus+=CalcdPhi(Phi, phi1);
1057           Nplus++;
1058
1059         }
1060
1061
1062       Phi+=((phiMinus+phiPlus)/(Nminus+Nplus+1));
1063
1064       if( TMath::Abs(CalcdPhi(phiDir, Phi)) > TMath::Pi()/2)
1065         Phi+=TMath::Pi();
1066
1067       Phi=CalcdPhi0To2pi(Phi, 0);
1068
1069       return kTRUE;
1070     }
1071
1072
1073   Double_t xmax = -100;
1074   Double_t xmin =  100;
1075
1076   xmean = x;
1077   ymean = y;
1078
1079
1080   for(Int_t ip=0; ip< Nev; ip++) {
1081     if(list.At(ip)<0) break;
1082
1083     TVector3* part = (TVector3*) GetAt(list.At(ip));
1084     if(!part) continue;
1085
1086     TVector3 vtmp(*part);
1087     GetLocalMom(vec, &vtmp);
1088
1089     Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
1090     Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
1091
1092     Double_t x1 = (xx-xmean)*TMath::Cos(Phi) +(yy-ymean)*TMath::Sin(Phi);
1093     if(x1 > xmax) xmax = x1;
1094     if(x1 < xmin) xmin = x1;
1095     //    printf("c3\n");
1096   }
1097
1098   if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
1099
1100
1101
1102
1103    while(Phi>TwoPi) Phi-=TwoPi;
1104
1105   return kTRUE;
1106 }
1107
1108 /*
1109 /////////////////////////////////////
1110
1111  TH2F *fhPhiEtaTrack;//!
1112
1113  TH1F *fhTMA_JAA[3];//!
1114  TH1F *fhTMA_JAp[3];//!
1115  TH1F *fhTMA_B1AA[3];//!
1116  TH1F *fhTMA_B1Ap[3];//!
1117  TH1F *fhTMA_B2AA[3];//!
1118  TH1F *fhTMA_B2Ap[3];//!
1119
1120
1121  Int_t     fPtJetNbin;
1122  TArrayD   fPtJetArray;
1123  Int_t     fPsiVsRNbin;
1124  Double_t  fRmax;
1125  Int_t     fPsiVsPhiNbin;
1126
1127  UInt_t   fFilterMask;
1128  Double_t fEtaTrackMax;
1129  Double_t fPtTrackMin;
1130  Double_t fPtTrackMax;
1131  Double_t fPtJmin;
1132  Double_t fPtJmax;
1133  Double_t fPtJ;
1134
1135  TVector3 fJvec;
1136 */
1137
1138 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AliAnalysisTaskJetShapeHM(TString comment, Bool_t kMC):TObject(),
1139 fComment(comment),
1140 fkMC(kMC),
1141 fkMCprod(0),
1142 fhEvent(0),
1143 fhMult(0),
1144 fhPtJ(0),
1145 fhPtJvsPtCorr(0),
1146 fhPsiVsR(0),
1147 fhPsiVsRPtJ(0),
1148 fhPhiEtaTrack(0),
1149 fhPsiVsR_MCtr(0),
1150 fhPsiVsRPtJ_MCtr(0),
1151 fhJetTrPtVsPartPt(0),
1152 fPtJetNbin(0),
1153 fPtJetArray(0),
1154 fPsiVsRNbin(0),
1155 fRmax(0),
1156 fPsiVsPhiNbin(0),
1157 fFilterMask(0),
1158 fEtaTrackMax(0),
1159 fPtTrackMin(0),
1160 fPtTrackMax(0),
1161 fPtJmin(0),
1162 fPtJmax(0),
1163 fPtJ(0),
1164 fJvec(0,0,0)
1165 {
1166
1167   for(Int_t i=0; i<3; i++)
1168     {
1169       fhTMA_JAA[i]=0;
1170       fhTMA_JAp[i]=0;
1171       fhTMA_B1AA[i]=0;
1172       fhTMA_B1Ap[i]=0;
1173       fhTMA_B2AA[i]=0;
1174       fhTMA_B2Ap[i]=0;
1175     }
1176
1177   for(Int_t i=0; i<3; i++)
1178     {
1179       for(Int_t j=0; j<2; j++)
1180         {
1181           fhPtresVsPt[i][j]=0;
1182           fhPhiresVsPhi[i][j]=0;
1183           fhEtaresVsEta[i][j]=0;
1184           fhRresVsPt[i][j]=0;
1185           fhDCAxy[i][j]=0;
1186           fhDCAz[i][j]=0;
1187           fhTrackPtEtaPhi[i][j]=0;
1188         }
1189     }
1190
1191
1192   fPtJetNbin = 0;
1193   fPtJetArray = 0;
1194
1195   SetRBins();
1196   SetPhiNbins();
1197   SetFilterMask();
1198   SetEtaTrackMax();
1199   SetPtTrackRange();
1200   SetPtJetRange();
1201
1202   fJvec.SetXYZ(0,0,0);
1203   fPtJ = 0;
1204 }
1205
1206
1207 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::~AliAnalysisTaskJetShapeHM()
1208 {
1209   /*
1210   if(fPtJetArray)
1211     delete [] fPtJetArray;
1212   fPtJetArray=0;
1213   */
1214 }
1215
1216
1217 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::SetPtJetBins(Int_t Nbin, Double_t *array)
1218 {
1219
1220   /*
1221   if(fPtJetArray)
1222     delete [] fPtJetArray;
1223   fPtJetArray=0;
1224
1225   fPtJetNbin = Nbin;
1226   fPtJetArray = new Double_t[fPtJetNbin+1];
1227   for(Int_t i=0; i<= fPtJetNbin; i++) fPtJetArray[i]= array[i];
1228   */
1229
1230
1231   fPtJetNbin = Nbin;
1232   fPtJetArray.Set(fPtJetNbin+1, array);
1233 }
1234
1235
1236 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::InitHistos()
1237 {
1238
1239   fhEvent    = Hist1D("hEvent" , 3  , -0.5,   2.5, "Event");
1240   fhMult     = Hist1D("hMult"  , 101, -0.5, 100.5, "N_{ch}");
1241
1242   fhPhiEtaTrack = Hist2D("hPhiEtaTrack", 100, 0, TMath::TwoPi(), 20, -1, 1., "#phi", "#eta");
1243
1244   if(fPtJetNbin<1) {
1245     Int_t Nbin = 13;
1246     Double_t array[]= {0., 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 130, 160, 200};
1247     SetPtJetBins(Nbin, array);
1248   }
1249
1250   fhPtJ      = Hist1D("hPtJ"  , fPtJetNbin, fPtJetArray.GetArray(), "Pt_{J} GeV/c");
1251   //  fhPsiVsR   = Hist1D("hPsiVsR", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1252   fhPsiVsR   = Hist3D("hPsiVsR", fPsiVsRNbin, 0., fRmax, 10, 0., 1., fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{Jt, frac}", "P_{J} GeV/c");
1253   fhPsiVsRPtJ   = Hist2D("hPsiVsRPtJ", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1254
1255   fhPtJvsPtCorr  = Hist2D("fhPtJvsPtCorr", fPtJetNbin, fPtJetArray.GetArray(), 100, -100, 100, "P_{tJ} GeV/c", "P_{tJ} - P_{tJ,corr} GeV/c" , 1);
1256
1257
1258
1259
1260   for(Int_t i=0; i<3; i++)
1261     {
1262       fhTMA_JAA[i]  = Hist1D(Form("fhTMA_JAA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1263       fhTMA_JAp[i]  = Hist1D(Form("fhTMA_JAp_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1264
1265       fhTMA_B1AA[i]  = Hist1D(Form("fhTMA_B1AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1266       fhTMA_B1Ap[i]  = Hist1D(Form("fhTMA_B1Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1267
1268       fhTMA_B2AA[i]  = Hist1D(Form("fhTMA_B2AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1269       fhTMA_B2Ap[i]  = Hist1D(Form("fhTMA_B2Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1270     }
1271
1272   if(fkMCprod)
1273     {
1274   fhPsiVsR_MCtr      = Hist1D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1275   fhPsiVsRPtJ_MCtr   = Hist2D("hPsiVsRPtJ_MCtr", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1276   fhJetTrPtVsPartPt  = Hist2D("fhJetTrPtVsPartPt", fPtJetNbin, fPtJetArray.GetArray(), 100, -1, 1, "P_{tJ,part} GeV/c", "1-P_{tJ,tr}/P_{tJ,part} GeV/c" , 1);
1277   const char *ch[2]={"m","p"};
1278       for(Int_t i=0; i<3; i++)
1279         {
1280           for(Int_t j=0; j<2; j++)
1281             {
1282           fhPtresVsPt[i][j]    = Hist2D(Form("hPtresVsPt%d%s",i,ch[j]), 100, 0, 100, 100, -0.5, 0.5, "P_{t}^{gen} GeV/c", "1-P_{t}^{rec}/P_{t}^{gen} GeV/c" );
1283           fhPhiresVsPhi[i][j]  = Hist2D(Form("hPhiresVsPhi%d%s",i,ch[j]), 600, 0, TMath::TwoPi(), 128, -0.2, 0.2, "#phi^{gen} rad", "#phi^{rec} - #phi^{gen} rad" );
1284           fhEtaresVsEta[i][j]  = Hist2D(Form("hEtaresVsEta%d%s",i,ch[j]), 200, -1, 1, 40, -0.2, 0.2, "#eta^{gen}", "#eta^{rec} - #eta^{gen}" );
1285           fhRresVsPt[i][j]     = Hist2D(Form("hRresVsPt%d%s",i,ch[j])  , 100, 0, 100, 500, -fRmax, fRmax, "P_{t, track}^{gen} GeV/c", "R^{gen}-R^{rec}" );
1286           fhDCAxy[i][j]        = Hist2D(Form("hDCAxy%d%s",i,ch[j]), 100, 0, 100, 300, -3, 3, Form("p_{t} of part. type %d%s [GeV/c]",i,ch[j]), "DCAxy [cm]");
1287           fhDCAz[i][j]         = Hist2D(Form("hDCAz%d%s",i,ch[j]) , 100, 0, 100, 300, -3, 3, Form("p_{t} of part. type %d%s [GeV/c]",i,ch[j]), "DCAz  [cm]") ;
1288           fhTrackPtEtaPhi[i][j]= Hist3D(Form("hTrackPtEtaPhi%d%s",i,ch[j]), 100, 0, 100, 100, -1, 1, 100, 0, TMath::TwoPi(),"P_{t} GeV/c ","#eta","#phi rad");
1289             }
1290         }
1291     }
1292
1293 }
1294
1295
1296
1297
1298
1299 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
1300 {
1301   list->Add(fhEvent);
1302   list->Add(fhMult);
1303   list->Add(fhPhiEtaTrack);
1304   list->Add(fhPtJ);
1305   list->Add(fhPsiVsR);
1306   list->Add(fhPsiVsRPtJ);
1307   list->Add(fhPtJvsPtCorr);
1308
1309
1310   for(Int_t i=0; i<3; i++)
1311     {
1312       list->Add(fhTMA_JAA[i]);
1313       list->Add(fhTMA_JAp[i]);
1314
1315       list->Add(fhTMA_B1AA[i]);
1316       list->Add(fhTMA_B1Ap[i]);
1317
1318       list->Add(fhTMA_B2AA[i]);
1319       list->Add(fhTMA_B2Ap[i]);
1320     }
1321
1322   if(fkMCprod)
1323     {
1324       list->Add(fhPsiVsR_MCtr);
1325       list->Add(fhPsiVsRPtJ_MCtr);
1326       list->Add(fhJetTrPtVsPartPt);
1327
1328       for(Int_t i=0; i<3; i++)
1329         {
1330           for(Int_t j=0; j<2; j++)
1331             {
1332               list->Add(fhPtresVsPt[i][j]);
1333               list->Add(fhPhiresVsPhi[i][j]);
1334               list->Add(fhEtaresVsEta[i][j]);
1335               list->Add(fhRresVsPt[i][j]);
1336               list->Add(fhDCAxy[i][j]);
1337               list->Add(fhDCAz[i][j]);
1338               list->Add(fhTrackPtEtaPhi[i][j]);
1339             }
1340         }
1341     }
1342 }
1343
1344
1345
1346
1347 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE,  AliAODJet *jet, Double_t DeltaPtJ)
1348 {
1349
1350
1351
1352   Int_t size = aodE->GetNumberOfTracks();
1353
1354   TClonesArray *arrP = 0;
1355
1356   if(fkMC || fkMCprod)
1357     {
1358       arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
1359      if(!arrP)
1360        {
1361          printf("ERROR: no Info about particles in AOD\n");
1362          return kFALSE;
1363        }
1364     }
1365
1366   if(fkMC)
1367      size = arrP->GetEntriesFast();
1368
1369   Int_t IndexArray[size];
1370   Int_t IndexArrayMC[size];
1371
1372   TClonesArray farray("TVector3", size);
1373
1374   Int_t counter = 0;
1375   Int_t counterMC = 0;
1376   Int_t counterAll = 0;
1377   Double_t pJ[3] = {0, 0, 0};
1378   Double_t pJmc[3] = {0, 0, 0};
1379
1380   AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1381   Double_t bfield = aodE->GetMagneticField();
1382   Double_t dca[2];
1383   Double_t cov[3];
1384
1385   TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
1386
1387
1388
1389
1390   for(int it = 0;it < size;++it){
1391     TVector3 vec;
1392     Int_t ch = -999;
1393     Int_t label = 0;
1394
1395
1396     if(fkMC)
1397       {
1398         AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
1399         if(!part)continue;
1400         if(!part->IsPhysicalPrimary())continue;
1401         if(part->Charge()==0)continue;
1402         vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1403       }
1404     else
1405       {
1406         AliAODTrack *tr = aodE->GetTrack(it);
1407         if(!tr) continue;
1408         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1409         label = tr->GetLabel();
1410         AliAODTrack tmp(*tr);
1411         tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1412         vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
1413         ch = tr->Charge();
1414       }
1415
1416
1417     if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1418     if(vec.Pt()< fPtTrackMin)continue;
1419     if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1420
1421
1422     new(farray[counterAll]) TVector3(vec);
1423     counterAll++;
1424
1425
1426     Double_t R = CalcR(vecJ, vec);
1427     if(R> fRmax) continue;
1428
1429     pJ[0]+=vec.Px();
1430     pJ[1]+=vec.Py();
1431     pJ[2]+=vec.Pz();
1432
1433     IndexArray[counter] = it;
1434     counter++;
1435
1436
1437
1438     // effics
1439     if(fkMCprod)
1440       {
1441         //        printf("A1\n");
1442         AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1443         if(!part)continue;
1444
1445         Int_t type = 0;
1446         if(!part->IsPhysicalPrimary()) type=1;
1447         if(label <0) type=2;
1448
1449         if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1450           //      printf("A4\n");
1451         if(ch==-1) ch=0;
1452
1453         Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1454         Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
1455         Double_t phiT = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1456
1457         fhPtresVsPt[type][ch]->Fill(part->Pt(), 1-vec.Pt()/part->Pt());
1458         fhPhiresVsPhi[type][ch]->Fill(phip, dphi);
1459         fhEtaresVsEta[type][ch]->Fill(part->Eta(), vec.Eta() - part->Eta());
1460         fhDCAxy[type][ch]->Fill(part->Pt(), dca[0]);
1461         fhDCAz[type][ch]->Fill( part->Pt(), dca[1]);
1462         fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), phiT);
1463
1464         TVector3 vecP(part->Px(), part->Py(), part->Pz());
1465         Double_t Rgen = CalcR(vecJ, vecP);
1466         fhRresVsPt[type][ch]->Fill(part->Pt(), Rgen - R);
1467
1468         pJmc[0]+=part->Px();
1469         pJmc[1]+=part->Py();
1470         pJmc[2]+=part->Pz();
1471
1472         IndexArrayMC[counterMC] = label;
1473         counterMC++;
1474           //      printf("A5\n");
1475       }
1476
1477   }
1478
1479   fhMult->Fill(counter);
1480   if(counter<1) return kFALSE;
1481
1482   fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1483
1484
1485   fPtJ =  TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1486   if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1487   fhPtJ->Fill(fPtJ);
1488
1489
1490
1491
1492
1493   fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
1494
1495   for(Int_t it = 0; it<counter; it++)
1496     {
1497       TVector3 vec;
1498
1499       if(fkMC)
1500         {
1501           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
1502           if(!part) continue;
1503           vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1504         }
1505       else
1506         {
1507           AliAODTrack *tr = aodE->GetTrack(IndexArray[it]);
1508           if(!tr) continue;
1509           AliAODTrack tmp(*tr);
1510           tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1511           if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
1512         }
1513
1514         Double_t R = CalcR(fJvec, vec);
1515         //      Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
1516         Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
1517         //      fhPsiVsR->Fill(R, probL);
1518         //      fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1519         fhPsiVsR->Fill(R,probL, fJvec.Mag());
1520         fhPsiVsRPtJ->Fill(R, fPtJ);
1521
1522         Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1523         fhPhiEtaTrack->Fill(phi, vec.Eta());
1524
1525     }
1526
1527
1528   if(fkMCprod)
1529     {
1530       TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
1531       Double_t fPtJMCtr=  fJvecMCtr.Perp();
1532
1533       fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
1534       for(Int_t it = 0; it<counterMC; it++)
1535         {
1536           TVector3 vec;
1537
1538           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(IndexArrayMC[it])));
1539           if(!part) continue;
1540           vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1541
1542           Double_t R = CalcR(fJvecMCtr, vec);
1543           //      Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1544           //      fhPsiVsR->Fill(R, probL);
1545           //      fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1546           //      Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1547           fhPsiVsR_MCtr->Fill(R);
1548           fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1549         }
1550     }
1551
1552   // ang. distr.
1553   //     fMyTool->Clean();
1554   AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1555
1556       fMyTool->SetNEntries(counterAll);
1557       fMyTool->SetVecJ(vecJ);
1558       fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1559       fMyTool->Init();
1560       Int_t scenario = 3;
1561
1562       /*
1563        to be added!!!!!!!!!
1564       fhTMA_B1AA[i]=0;
1565       fhTMA_B2AA[i]=0;
1566       */
1567
1568       for(Int_t l=0; l<3; l++)
1569         {
1570           Double_t phiA, phi1;
1571
1572           //      Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
1573           //      Double_t ptRJ  = fMyTool->GetPRecJ(l,0).Pt();
1574           Int_t N1 = fMyTool->GetSizeJ(l,1);
1575           Double_t dphi = -999;
1576
1577
1578           if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0)  , vecJ, phiA, scenario))
1579             continue;
1580
1581           //      f2JEvent->SetPhiJL(l,0, phiA);
1582
1583           if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1)  , vecJ, phi1, scenario))
1584             {
1585               fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1586             }
1587
1588
1589           for(Int_t ip =0; ip<N1; ip++)
1590             {
1591               phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1592               dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1593               fhTMA_JAp[l]->Fill(dphi);
1594             }
1595
1596
1597           Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1598           for(Int_t ip =0; ip<NB1; ip++)
1599             {
1600               phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1601               dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1602               fhTMA_B1Ap[l]->Fill(dphi);
1603             }
1604
1605           Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1606           for(Int_t ip =0; ip<NB2; ip++)
1607             {
1608               phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1609               dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1610               fhTMA_B2Ap[l]->Fill(dphi);
1611             }
1612
1613         }
1614
1615     fhEvent->Fill(1);
1616
1617     delete fMyTool;
1618
1619   return kTRUE;
1620 }
1621
1622
1623 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
1624 {
1625
1626   Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
1627   Double_t deta = v1.Eta() - v2.Eta();
1628   Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
1629   return RB;
1630 }
1631
1632
1633 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1634 {
1635
1636   while(phi1<0) phi1+=TMath::TwoPi();
1637   while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1638
1639   while(phi2<0) phi2+=TMath::TwoPi();
1640   while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
1641
1642   Double_t dphi = phi1- phi2;
1643   if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
1644   if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
1645
1646   return dphi;
1647 }
1648
1649
1650
1651
1652
1653
1654
1655 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax,  const char* xLabel, Int_t color, const char* yLabel)
1656 {
1657 // create a 1D histogram
1658
1659   TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
1660   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1661   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1662   res->SetLineColor(color);
1663   return res;
1664 }
1665
1666
1667 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray,  const char* xLabel, Int_t color, const char* yLabel)
1668 {
1669 // create a 1D histogram
1670
1671   TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
1672   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1673   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1674   res->SetLineColor(color);
1675   return res;
1676 }
1677
1678 TH2F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel, const char* yLabel, Int_t color)
1679 {
1680 // create a 2D histogram
1681
1682   TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
1683   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1684   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1685   //  res->SetMarkerStyle(kFullCircle);
1686   //  res->SetOption("E");
1687   res->SetLineColor(color);
1688   //  fOutputList->Add(res);
1689   return res;
1690 }
1691
1692
1693 TH2F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t *yArray, const char* xLabel, const char* yLabel, Int_t color, const char* zLabel)
1694 {
1695 // create a 2D histogram
1696
1697   TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
1698   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1699   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1700   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1701   //  res->SetMarkerStyle(kFullCircle);
1702   //  res->SetOption("E");
1703   res->SetLineColor(color);
1704   //  fOutputList->Add(res);
1705   return res;
1706 }
1707
1708 TH2F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist2D(const char* name, Int_t nBinsx, Double_t *yArrax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel, const char* yLabel, Int_t color, const char* zLabel)
1709 {
1710 // create a 2D histogram
1711
1712   TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
1713   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1714   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1715   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1716   //  res->SetMarkerStyle(kFullCircle);
1717   //  res->SetOption("E");
1718   res->SetLineColor(color);
1719   //  fOutputList->Add(res);
1720   return res;
1721 }
1722
1723 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1724                                                   Int_t nBinsy, Double_t yMin, Double_t yMax,
1725                                                   Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1726 {
1727 // create a 3D histogram
1728
1729   TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
1730   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1731   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1732   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1733   //  res->SetMarkerStyle(kFullCircle);
1734   //  res->SetOption("E");
1735   res->SetLineColor(color);
1736   //  fOutputList->Add(res);
1737   return res;
1738 }
1739
1740 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1741                                                   Int_t nBinsy, Double_t yMin, Double_t yMax,
1742                                                   Int_t nBinsz, Double_t *zArr, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1743 {
1744 // create a 3D histogram
1745
1746   Double_t xArr[nBinsx+1], yArr[nBinsy+1];
1747
1748   for(Int_t i=0; i<=nBinsx; i++) xArr[i]=xMin+i*(xMax-xMin)/nBinsx;
1749   for(Int_t i=0; i<=nBinsy; i++) yArr[i]=yMin+i*(yMax-yMin)/nBinsy;
1750
1751   TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xArr, nBinsy, yArr, nBinsz, zArr);
1752   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1753   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1754   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1755   //  res->SetMarkerStyle(kFullCircle);
1756   //  res->SetOption("E");
1757   res->SetLineColor(color);
1758   //  fOutputList->Add(res);
1759   return res;
1760 }
1761
1762
1763 //////////////////////////////////////
1764
1765
1766 /*
1767 //________________________________________________________________________
1768 void AnalysisJetMain::ConnectInputData(Option_t *)
1769 {
1770   // Connect ESD
1771   // Called once
1772
1773   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1774   //   if (!fESD) {
1775   //     AliError("ESD not available");
1776   //     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1777   //     fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1778
1779   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1780   //  assume that the AOD is in the general output...
1781   fAODOut  = AODEvent();
1782
1783
1784 }
1785 */
1786
1787 void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1788 {
1789    fJetBranchName[0] = branch1;
1790    fJetBranchName[1] = branch2;
1791 }
1792
1793 //________________________________________________________________________
1794 void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1795 {
1796   //printf("Open1\n");
1797   //       const char *nameF = OpenFile(1)->GetName();
1798   //printf("Open2 %s\n",nameF);
1799
1800   //  fTriggerAnalysis = new AliTriggerAnalysis();
1801   //  if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1802
1803
1804   fOutputList = new TList();
1805   fOutputList->SetOwner();
1806
1807   fhPtJL  = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1808   fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1809
1810
1811
1812   printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1813
1814
1815   fanalJetSubStrHM->SetFilterMask(fFilterMask);
1816   if(fkMC)   fanalJetSubStrHM->MCprod();
1817   fanalJetSubStrHM->InitHistos();
1818   fanalJetSubStrHM->AddToList(fOutputList);
1819
1820
1821   if(fkMC)
1822     {
1823
1824       printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1825       fanalJetSubStrMCHM->InitHistos();
1826       fanalJetSubStrMCHM->AddToList(fOutputList);
1827
1828       for(Int_t i=0; i<3; i++)
1829         {
1830           fhPtresVsPt[i]    = Hist2D(Form("hPtresVsPt%d",i), 100, 0, 100, 100, -0.5, 0.5, "P_{t}^{gen} GeV/c", "1-P_{t}^{rec}/P_{t}^{gen} GeV/c" );      fOutputList->Add(fhPtresVsPt[i]);
1831           fhPhiresVsPhi[i]  = Hist2D(Form("hPhiresVsPhi%d",i), 600, 0, TMath::TwoPi(), 128, -0.2, 0.2, "#phi^{gen} rad", "#phi^{rec} - #phi^{gen} rad" );      fOutputList->Add(fhPhiresVsPhi[i]);
1832           fhEtaresVsEta[i]  = Hist2D(Form("hEtaresVsEta%d",i), 200, -1, 1, 40, -0.2, 0.2, "#eta^{gen}", "#eta^{rec} - #eta^{gen}" );      fOutputList->Add(fhEtaresVsEta[i]);
1833           fhDCAxy[i]        = Hist2D(Form("hDCAxy%d",i), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]"); fOutputList->Add(fhDCAxy[i]);
1834           fhDCAz[i]         = Hist2D(Form("hDCAz%d",i) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ; fOutputList->Add(fhDCAz[i]);
1835           fhTrackPtEtaPhi[i]= Hist3D(Form("hTrackPtEtaPhi%d",i), 100, 0, 100, 100, -1, 1, 100, 0, TMath::TwoPi(),"P_{t} GeV/c ","#eta","#phi rad");  fOutputList->Add(fhTrackPtEtaPhi[i]);
1836         }
1837     }
1838
1839
1840
1841
1842    printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1843
1844
1845
1846   PostData(1, fOutputList);
1847
1848   /*
1849    OpenFile(1);
1850
1851   fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1852   fOutputTree->Branch("event",&f2JEvent);
1853   */
1854
1855 //    PostData(1, fOutputTree);
1856 }
1857
1858 /*
1859 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1860 {
1861
1862   // Fetch the aod also from the input in,
1863   // have todo it in notify
1864
1865
1866   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1867   //  assume that the AOD is in the general output...
1868   fAODOut  = AODEvent();
1869
1870   if(fNonStdFile.Length()!=0){
1871     // case that we have an AOD extension we need can fetch the jets from the extended output
1872     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1873     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1874     if(!fAODExtension){
1875       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1876     }
1877   }
1878
1879
1880
1881   return kTRUE;
1882 }
1883 */
1884
1885
1886
1887
1888 //________________________________________________________________________
1889 void AliAnalysisTaskJetShape::UserExec(Option_t *)
1890 {
1891   //   return;
1892 //    if(f2JEvent)
1893 //      delete f2JEvent;
1894 //    f2JEvent = new Ali2JEvent();
1895 //    return;
1896
1897   //  return;
1898  if(fDebug)  Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
1899   if(fDebug) printf("NEW EVENT 0\n");
1900
1901   if(!IsGoodEvent()) return;
1902
1903
1904   //   printf("\n\n\n NEW EVENT\n");
1905
1906    AliAODEvent* aodE = 0;
1907    if(!aodE){
1908      if(!fESD)aodE = fAODIn;
1909      else aodE = fAODOut;}
1910
1911    /*
1912    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1913 if(!aodH){
1914     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1915     return;
1916   }
1917    */
1918    if(fDebug) {
1919      printf("NEW EVENT 1:  Number of tracks = %d\n",aodE->GetNumberOfTracks());
1920      //     aodE->GetList()->Print();
1921      //     printf("Look for %s\n",fJetBranchName[0].Data());
1922    }
1923
1924
1925   // centrality selection
1926    AliCentrality *cent = 0x0;
1927    Double_t centrality = 0.;
1928
1929    if(fESD) {cent = fESD->GetCentrality();
1930      if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1931    else     centrality=aodE->GetHeader()->GetCentrality();
1932
1933
1934    if(!fkIsPbPb) {
1935      centrality=1.;
1936    }
1937
1938 //  if(fDebug) printf("NEW EVENT 2\n");
1939
1940    Bool_t  fUseAOD049 = kFALSE;// kTRUE;//
1941       if(fUseAOD049&&centrality>=0){
1942         Float_t v0=0;
1943         AliAODVZERO* aodV0 = aodE->GetVZEROData();
1944         v0+=aodV0->GetMTotV0A();
1945         v0+=aodV0->GetMTotV0C();
1946         if(centrality==0 && v0 < 19500) return ;//filtering issue
1947         Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
1948         Float_t val= 1.30552 +  0.147931 * v0;
1949         Float_t tklSigma[101]={176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86, 120.788, 115.611, 113.172, 110.496,
1950             109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654, 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1951             68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224, 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114,
1952             43.2083, 41.3065, 40.1863, 38.5255, 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398, 26.6488, 25.0183,
1953             25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235, 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173,
1954             12.9504, 12.9504, 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544,
1955             13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1956         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] )
1957           return; //outlier
1958       }
1959
1960
1961     if(fDebug)
1962       {
1963         printf("centrality: %f\n", centrality);
1964         aodE->Print();
1965         aodE->GetList()->Print();
1966       }
1967
1968
1969     if (centrality < fCentMin || centrality > fCentMax){
1970     //    PostData(1, fOutputList);
1971     return;
1972     }
1973
1974    //   fhNEvent->Fill(0);
1975
1976    // accepted events
1977    // -- end event selection --
1978
1979
1980    // get background
1981    AliAODJetEventBackground* externalBackground = 0;
1982    AliAODJetEventBackground* externalBackgroundMC = 0;
1983
1984
1985    Float_t rho = 0;
1986    Float_t rho_MC=0;
1987
1988    if(fkIsPbPb)
1989      {
1990        if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
1991          externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
1992          if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1993        }
1994        if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
1995          externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
1996          if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1997        }
1998        if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
1999          externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
2000          if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2001        }
2002        if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2003          externalBackgroundMC =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
2004          if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2005        }
2006        if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2007          externalBackgroundMC =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
2008          if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2009        }
2010        if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2011          externalBackgroundMC =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
2012          if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2013        }
2014        if(externalBackground)rho = externalBackground->GetBackground(0);
2015        if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
2016      }
2017
2018
2019
2020   if(fkMC)
2021     {
2022       AliAODVertex* primVtx = aodE->GetPrimaryVertex();
2023       Double_t bfield = aodE->GetMagneticField();
2024
2025       TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2026       if(!arrP)
2027         {
2028           printf("ERROR: no Info about particles in AOD\n");
2029           return;
2030         }
2031
2032       for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2033         {
2034           AliAODTrack *tr = aodE->GetTrack(it);
2035           if(!tr) continue;
2036           if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2037           if(TMath::Abs(tr->Eta())>1.) continue;
2038           Int_t label = TMath::Abs(tr->GetLabel());
2039
2040           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2041           if(!part)continue;
2042           //      if(!part->IsPhysicalPrimary())continue;
2043           //      if(part->Charge()==0)continue;
2044           //      vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2045
2046
2047
2048           Double_t dca[2];
2049           Double_t cov[3];
2050
2051           AliAODTrack tmp(*tr);
2052           tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2053
2054           Int_t type = 0;
2055           if(!part->IsPhysicalPrimary()) type=1;
2056           if(label<0) type =2;
2057
2058           fhPtresVsPt[type]->Fill(part->Pt(), 1-tr->Pt()/part->Pt());
2059           fhPhiresVsPhi[type]->Fill(part->Phi(), tr->Phi() - part->Phi());
2060           fhEtaresVsEta[type]->Fill(part->Eta(), tr->Eta() - part->Eta());
2061           fhDCAxy[type]->Fill(part->Pt(), dca[0]);
2062           fhDCAz[type]->Fill( part->Pt(), dca[1]);
2063           fhTrackPtEtaPhi[type]->Fill(tr->Pt(), tr->Eta(), tr->Phi());
2064         }
2065
2066
2067     }
2068
2069
2070
2071
2072
2073    //   printf("rho = %f\n",rho);
2074
2075     //   return;
2076
2077    TClonesArray *Jets[2];
2078    Jets[0]=0;
2079    Jets[1]=0;
2080    if(fAODOut&&!Jets[0]){
2081    Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
2082    Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
2083    if(fAODExtension && !Jets[0]){
2084    Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
2085    Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
2086    if(fAODIn&&!Jets[0]){
2087    Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
2088    Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
2089
2090
2091    if(!Jets[0]) {
2092      Printf("no friend!!!\n");
2093      return;
2094    }
2095   Int_t nJ = Jets[0]->GetEntries();
2096   Float_t ptmax = 0.;
2097   Int_t   imax  = -1;
2098
2099   if(fDebug) {
2100     printf("NEW EVENT 3:  Number of Rec. jets  %d\n",nJ);
2101   }
2102
2103 //
2104 // Find highest pT jet with pt > 20 GeV
2105 //
2106 //  fPtJcorrMin=0;
2107
2108   Float_t etaJmax = 0.4;
2109
2110      for (Int_t i = 0; i < nJ; i++) {
2111        AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2112        if (!jet) continue;
2113        //              jet->Print("0");
2114        Float_t ptJ  = jet->Pt();
2115        Float_t etaJ = TMath::Abs(jet->Eta());
2116
2117        Float_t area = jet->EffectiveAreaCharged();
2118        Float_t ptJcorr=ptJ-rho*area;
2119
2120        if ((ptJcorr > fPtJcorrMin) && (ptJcorr  > ptmax) && etaJ < etaJmax) {
2121          ptmax = ptJcorr;
2122          imax = i;
2123        }
2124
2125      }
2126
2127
2128   if(fDebug) {
2129     printf("NEW EVENT 4:\n");
2130   }
2131
2132      TVector3 vecJ;
2133
2134      AliAODJet* jetL = 0;
2135
2136      Float_t areaJL, ptJLcorr;
2137
2138      if (imax != -1)  {
2139
2140        jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2141        if(jetL){
2142
2143          areaJL   = jetL->EffectiveAreaCharged();
2144          ptJLcorr  = jetL->Pt()-rho*areaJL;
2145
2146          fhPtJL->Fill(ptJLcorr);
2147          fhAreaJL->Fill(areaJL);
2148          vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2149          fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
2150
2151          if(fDebug) {
2152            Printf("Leading Jet");
2153            jetL->Print("0");
2154          }
2155
2156        }
2157      }
2158
2159
2160      if(!fkMC)
2161        {
2162          PostData(1, fOutputList);
2163          if(fDebug)  Printf("End of event #%5d", (Int_t) fEntry);
2164          return;
2165        }
2166
2167
2168
2169
2170      nJ = Jets[1]->GetEntries();
2171   if(fDebug) {
2172     printf("NEW EVENT 5:  Number of Rec. jets  %d\n",nJ);
2173   }
2174
2175      ptmax = 0;
2176      imax = -1;
2177      Float_t areaJL_MC=0;
2178
2179      for (Int_t i = 0; i < nJ; i++) {
2180        AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2181        if (!jet) continue;
2182        Float_t ptJ1  = jet->Pt();
2183        Float_t etaJ1 = TMath::Abs(jet->Eta());
2184
2185        areaJL_MC = jet->EffectiveAreaCharged();
2186        Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2187
2188        //    jet->Print();
2189
2190        if ((ptJcorr > fPtJcorrMin) && (ptJcorr  > ptmax) && etaJ1 < etaJmax) {
2191          ptmax = ptJcorr;
2192          imax = i;
2193        }
2194
2195      }
2196
2197   if(fDebug) {
2198     printf("NEW EVENT 6:\n");
2199   }
2200
2201
2202      AliAODJet* jetMCL=0;
2203
2204     if (imax != -1)  {
2205       jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
2206       if(jetMCL){
2207         fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2208       }
2209     }
2210
2211
2212
2213
2214
2215      PostData(1, fOutputList);
2216
2217     if(fDebug)  Printf("End of event #%5d", (Int_t) fEntry);
2218
2219      return;
2220
2221
2222
2223
2224 }
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238 Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2239 {
2240   while(phi1<0) phi1+=TMath::TwoPi();
2241   while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2242
2243   while(phi2<0) phi2+=TMath::TwoPi();
2244   while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
2245
2246   Double_t dphi = phi1- phi2;
2247   if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
2248   if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
2249
2250   //  Double_t dphi = phi1- phi2;
2251   //  while(dphi<0) dphi+=TMath::TwoPi();
2252   //  while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2253   return dphi;
2254 }
2255
2256 Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2257 {
2258
2259   Double_t dphi = CalcdPhi(phi1, phi2);
2260   while(dphi<0) dphi+=TMath::TwoPi();
2261   while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2262   return dphi;
2263 }
2264
2265
2266 Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2267 {
2268
2269   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2270   //   if (!fESD) {
2271   //     AliError("ESD not available");
2272   //     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2273   //     fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2274
2275   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2276   //  assume that the AOD is in the general output...
2277   fAODOut  = AODEvent();
2278
2279
2280    static AliAODEvent* aod = 0;
2281    // take all other information from the aod we take the tracks from
2282    if(!aod){
2283      if(!fESD)aod = fAODIn;
2284      else aod = fAODOut;}
2285
2286
2287
2288    if(fNonStdFile.Length()!=0){
2289     // case that we have an AOD extension we need can fetch the jets from the extended output
2290      AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2291      fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2292    }
2293
2294    if(fDebug){
2295     if(fAODIn) printf("fAODIn\n");
2296     if(fAODOut) printf("fAODOut\n");
2297     if(fAODExtension) printf("fAODExtension\n");
2298    }
2299
2300
2301    /*
2302    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2303 if(!aodH){
2304     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2305       return kFALSE;
2306   }
2307    */
2308
2309    // -- event selection --
2310
2311    // physics selection
2312    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2313    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2314    //     cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
2315    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2316       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2317       return kFALSE;
2318    }
2319
2320    // vertex selection
2321    //   if(!aod){
2322    //     if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2323    //     fhNEvent->Fill(3);
2324    //      PostData(1, fOutputList);
2325    //     return kFALSE;
2326    //   }
2327
2328    AliAODVertex* primVtx = aod->GetPrimaryVertex();
2329
2330    if(!primVtx){
2331      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
2332      return kFALSE;
2333    }
2334
2335    Int_t nTracksPrim = primVtx->GetNContributors();
2336    if ((nTracksPrim < fVtxMinContrib) ||
2337          (primVtx->GetZ() < fVtxZMin) ||
2338          (primVtx->GetZ() > fVtxZMax) ){
2339       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2340       return  kFALSE;
2341    }
2342
2343    /*
2344    // event class selection (from jet helper task)
2345    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
2346    if(fDebug) Printf("Event class %d", eventClass);
2347    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
2348       return  kFALSE;
2349    }
2350    */
2351     return kTRUE;
2352 }
2353
2354
2355
2356
2357
2358
2359 //________________________________________________________________________
2360 void AliAnalysisTaskJetShape::Terminate(Option_t *)
2361 {
2362   // Draw result to the screen
2363   // Called once at the end of the query
2364  //    fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2365
2366   printf("MyTaskTestTrigger::Terminate()\n\n\n");
2367
2368
2369   //  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2370   //  fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2371
2372 }
2373
2374
2375
2376 /*
2377 Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd,  Double_t Vxyz[3], Int_t type)
2378 {
2379   // type
2380   //     0 -> SPD vtx
2381   //     1 -> ESD vtx
2382   const AliESDVertex* vtx = 0;
2383
2384   if(type==0) {
2385     vtx = esd->GetPrimaryVertexSPD();
2386     if(!vtx) return kFALSE;
2387     if(vtx->GetNContributors() < 1) return kFALSE;
2388     TString vtxTyp = vtx->GetTitle();
2389     if (vtxTyp.Contains("vertexer: Z")) {
2390       if (vtx->GetDispersion()>0.04) return kFALSE;
2391       if (vtx->GetZRes()>0.25) return kFALSE;
2392     }
2393
2394   }
2395   if(type==1) {
2396     vtx = esd->GetPrimaryVertexTracks();
2397     if(!vtx) return kFALSE;
2398     if(vtx->GetNContributors() < 1) return kFALSE;
2399   }
2400
2401
2402    Vxyz[0] = vtx->GetXv();
2403    Vxyz[1] = vtx->GetYv();
2404    Vxyz[2] = vtx->GetZv();
2405
2406    return kTRUE;
2407 }
2408 */
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418 TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax,  const char* xLabel, Int_t color)
2419 {
2420 // create a 1D histogram
2421
2422   TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
2423   if (xLabel) h->GetXaxis()->SetTitle(xLabel);
2424   //  res->SetMarkerStyle(kFullCircle);
2425   //  res->SetOption("E");
2426   h->SetLineColor(color);
2427   //  fOutputList->Add(h);
2428   return h;
2429 }
2430
2431
2432 TH2F *AliAnalysisTaskJetShape::Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel, const char* yLabel, Int_t color)
2433 {
2434 // create a 2D histogram
2435
2436   TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
2437   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2438   if (xLabel) res->GetYaxis()->SetTitle(yLabel);
2439   //  res->SetMarkerStyle(kFullCircle);
2440   //  res->SetOption("E");
2441   res->SetLineColor(color);
2442   //  fOutputList->Add(res);
2443   return res;
2444 }
2445
2446
2447 TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
2448                                                   Int_t nBinsy, Double_t yMin, Double_t yMax,
2449                                                   Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
2450 {
2451 // create a 3D histogram
2452
2453   TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
2454   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2455   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
2456   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
2457   //  res->SetMarkerStyle(kFullCircle);
2458   //  res->SetOption("E");
2459   res->SetLineColor(color);
2460   //  fOutputList->Add(res);
2461   return res;
2462 }
2463
2464
2465
2466
2467 //#endif
2468
2469
2470
2471
2472
2473
2474
2475