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, 100, 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      = Hist3D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, 100, 0., 1., fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{Jt, frac}", "P_{J} GeV/c");
1275   //  fhPsiVsR_MCtr      = Hist1D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1276   fhPsiVsRPtJ_MCtr   = Hist2D("hPsiVsRPtJ_MCtr", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1277   fhJetTrPtVsPartPt  = Hist2D("fhJetTrPtVsPartPt", fPtJetNbin, fPtJetArray.GetArray(), 100, -1, 1, "P_{tJ,part} GeV/c", "1-P_{tJ,tr}/P_{tJ,part} GeV/c" , 1);
1278   const char *ch[2]={"m","p"};
1279       for(Int_t i=0; i<3; i++)
1280         {
1281           for(Int_t j=0; j<2; j++)
1282             {
1283           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" );
1284           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" );
1285           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}" );
1286           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}" );
1287           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]");
1288           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]") ;
1289           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");
1290             }
1291         }
1292     }
1293
1294 }
1295
1296
1297
1298
1299
1300 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
1301 {
1302   list->Add(fhEvent);
1303   list->Add(fhMult);
1304   list->Add(fhPhiEtaTrack);
1305   list->Add(fhPtJ);
1306   list->Add(fhPsiVsR);
1307   list->Add(fhPsiVsRPtJ);
1308   list->Add(fhPtJvsPtCorr);
1309
1310
1311   for(Int_t i=0; i<3; i++)
1312     {
1313       list->Add(fhTMA_JAA[i]);
1314       list->Add(fhTMA_JAp[i]);
1315
1316       list->Add(fhTMA_B1AA[i]);
1317       list->Add(fhTMA_B1Ap[i]);
1318
1319       list->Add(fhTMA_B2AA[i]);
1320       list->Add(fhTMA_B2Ap[i]);
1321     }
1322
1323   if(fkMCprod)
1324     {
1325       list->Add(fhPsiVsR_MCtr);
1326       list->Add(fhPsiVsRPtJ_MCtr);
1327       list->Add(fhJetTrPtVsPartPt);
1328
1329       for(Int_t i=0; i<3; i++)
1330         {
1331           for(Int_t j=0; j<2; j++)
1332             {
1333               list->Add(fhPtresVsPt[i][j]);
1334               list->Add(fhPhiresVsPhi[i][j]);
1335               list->Add(fhEtaresVsEta[i][j]);
1336               list->Add(fhRresVsPt[i][j]);
1337               list->Add(fhDCAxy[i][j]);
1338               list->Add(fhDCAz[i][j]);
1339               list->Add(fhTrackPtEtaPhi[i][j]);
1340             }
1341         }
1342     }
1343 }
1344
1345
1346
1347
1348 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE,  AliAODJet *jet, Double_t DeltaPtJ)
1349 {
1350
1351
1352
1353   Int_t size = aodE->GetNumberOfTracks();
1354
1355   TClonesArray *arrP = 0;
1356
1357   if(fkMC || fkMCprod)
1358     {
1359       arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
1360      if(!arrP)
1361        {
1362          printf("ERROR: no Info about particles in AOD\n");
1363          return kFALSE;
1364        }
1365     }
1366
1367   if(fkMC)
1368      size = arrP->GetEntriesFast();
1369
1370   Int_t IndexArray[size];
1371   Int_t IndexArrayMC[size];
1372
1373   TClonesArray farray("TVector3", size);
1374
1375   Int_t counter = 0;
1376   Int_t counterMC = 0;
1377   Int_t counterAll = 0;
1378   Double_t pJ[3] = {0, 0, 0};
1379   Double_t pJmc[3] = {0, 0, 0};
1380
1381   AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1382   Double_t bfield = aodE->GetMagneticField();
1383   Double_t dca[2] = {0., 0.};
1384   Double_t cov[3] = {0., 0., 0.};
1385
1386   TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
1387
1388
1389
1390
1391   for(int it = 0;it < size;++it){
1392     TVector3 vec;
1393     Int_t ch = -999;
1394     Int_t label = 0;
1395
1396
1397     if(fkMC)
1398       {
1399         AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
1400         if(!part)continue;
1401         if(!part->IsPhysicalPrimary())continue;
1402         if(part->Charge()==0)continue;
1403         vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1404       }
1405     else
1406       {
1407         AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodE->GetTrack(it));
1408         if(!tr) AliFatal("Not a standard AOD");
1409         if(!tr) continue;
1410         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1411         label = tr->GetLabel();
1412         AliAODTrack tmp(*tr);
1413         tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1414         vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
1415         ch = tr->Charge();
1416       }
1417
1418
1419     if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1420     if(vec.Pt()< fPtTrackMin)continue;
1421     if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1422
1423
1424     new(farray[counterAll]) TVector3(vec);
1425     counterAll++;
1426
1427
1428     Double_t R = CalcR(vecJ, vec);
1429     if(R> fRmax) continue;
1430
1431     pJ[0]+=vec.Px();
1432     pJ[1]+=vec.Py();
1433     pJ[2]+=vec.Pz();
1434
1435     IndexArray[counter] = it;
1436     counter++;
1437
1438
1439
1440     // effics
1441     if(fkMCprod)
1442       {
1443         //        printf("A1\n");
1444         AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1445         if(!part)continue;
1446
1447         Int_t type = 0;
1448         if(!part->IsPhysicalPrimary()) type=1;
1449         if(label <0) type=2;
1450
1451         if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1452           //      printf("A4\n");
1453         if(ch==-1) ch=0;
1454
1455         Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1456         Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
1457         Double_t phiT = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1458
1459         fhPtresVsPt[type][ch]->Fill(part->Pt(), 1-vec.Pt()/part->Pt());
1460         fhPhiresVsPhi[type][ch]->Fill(phip, dphi);
1461         fhEtaresVsEta[type][ch]->Fill(part->Eta(), vec.Eta() - part->Eta());
1462         fhDCAxy[type][ch]->Fill(part->Pt(), dca[0]);
1463         fhDCAz[type][ch]->Fill( part->Pt(), dca[1]);
1464         fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), phiT);
1465
1466         TVector3 vecP(part->Px(), part->Py(), part->Pz());
1467         Double_t Rgen = CalcR(vecJ, vecP);
1468         fhRresVsPt[type][ch]->Fill(part->Pt(), Rgen - R);
1469
1470         pJmc[0]+=part->Px();
1471         pJmc[1]+=part->Py();
1472         pJmc[2]+=part->Pz();
1473
1474         IndexArrayMC[counterMC] = label;
1475         counterMC++;
1476           //      printf("A5\n");
1477       }
1478
1479   }
1480
1481   fhMult->Fill(counter);
1482   if(counter<1) return kFALSE;
1483
1484   fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1485
1486
1487   fPtJ =  TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1488   if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1489   fhPtJ->Fill(fPtJ);
1490
1491
1492
1493
1494
1495   fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
1496
1497   for(Int_t it = 0; it<counter; it++)
1498     {
1499       TVector3 vec;
1500
1501       if(fkMC)
1502         {
1503           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
1504           if(!part) continue;
1505           vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1506         }
1507       else
1508         {
1509           AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodE->GetTrack(IndexArray[it]));
1510           if(!tr) AliFatal("Not a standard AOD");
1511           if(!tr) continue;
1512           AliAODTrack tmp(*tr);
1513           tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1514           if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
1515         }
1516
1517         Double_t R = CalcR(fJvec, vec);
1518         //      Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
1519         Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
1520         //      fhPsiVsR->Fill(R, probL);
1521         //      fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1522         fhPsiVsR->Fill(R,probL, fJvec.Mag());
1523         fhPsiVsRPtJ->Fill(R, fPtJ);
1524
1525         Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1526         fhPhiEtaTrack->Fill(phi, vec.Eta());
1527
1528     }
1529
1530
1531   if(fkMCprod)
1532     {
1533       TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
1534       Double_t fPtJMCtr=  fJvecMCtr.Perp();
1535
1536       fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
1537       for(Int_t it = 0; it<counterMC; it++)
1538         {
1539           TVector3 vec;
1540
1541           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(IndexArrayMC[it])));
1542           if(!part) continue;
1543           vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1544
1545           Double_t R = CalcR(fJvecMCtr, vec);
1546           Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1547           //      fhPsiVsR->Fill(R, probL);
1548           //      fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1549           //Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1550           fhPsiVsR_MCtr->Fill(R,probL, fJvecMCtr.Mag());
1551           fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1552         }
1553     }
1554
1555   // ang. distr.
1556   //     fMyTool->Clean();
1557   AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1558
1559       fMyTool->SetNEntries(counterAll);
1560       fMyTool->SetVecJ(vecJ);
1561       fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1562       fMyTool->Init();
1563       Int_t scenario = 3;
1564
1565       /*
1566        to be added!!!!!!!!!
1567       fhTMA_B1AA[i]=0;
1568       fhTMA_B2AA[i]=0;
1569       */
1570
1571       for(Int_t l=0; l<3; l++)
1572         {
1573           Double_t phiA, phi1;
1574
1575           //      Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
1576           //      Double_t ptRJ  = fMyTool->GetPRecJ(l,0).Pt();
1577           Int_t N1 = fMyTool->GetSizeJ(l,1);
1578           Double_t dphi = -999;
1579
1580
1581           if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0)  , vecJ, phiA, scenario))
1582             continue;
1583
1584           //      f2JEvent->SetPhiJL(l,0, phiA);
1585
1586           if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1)  , vecJ, phi1, scenario))
1587             {
1588               fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1589             }
1590
1591
1592           for(Int_t ip =0; ip<N1; ip++)
1593             {
1594               phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1595               dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1596               fhTMA_JAp[l]->Fill(dphi);
1597             }
1598
1599
1600           Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1601           for(Int_t ip =0; ip<NB1; ip++)
1602             {
1603               phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1604               dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1605               fhTMA_B1Ap[l]->Fill(dphi);
1606             }
1607
1608           Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1609           for(Int_t ip =0; ip<NB2; ip++)
1610             {
1611               phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1612               dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1613               fhTMA_B2Ap[l]->Fill(dphi);
1614             }
1615
1616         }
1617
1618     fhEvent->Fill(1);
1619
1620     delete fMyTool;
1621
1622   return kTRUE;
1623 }
1624
1625
1626 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
1627 {
1628
1629   Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
1630   Double_t deta = v1.Eta() - v2.Eta();
1631   Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
1632   return RB;
1633 }
1634
1635
1636 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1637 {
1638
1639   while(phi1<0) phi1+=TMath::TwoPi();
1640   while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1641
1642   while(phi2<0) phi2+=TMath::TwoPi();
1643   while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
1644
1645   Double_t dphi = phi1- phi2;
1646   if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
1647   if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
1648
1649   return dphi;
1650 }
1651
1652
1653
1654
1655
1656
1657
1658 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)
1659 {
1660 // create a 1D histogram
1661
1662   TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
1663   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1664   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1665   res->SetLineColor(color);
1666   return res;
1667 }
1668
1669
1670 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray,  const char* xLabel, Int_t color, const char* yLabel)
1671 {
1672 // create a 1D histogram
1673
1674   TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
1675   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1676   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1677   res->SetLineColor(color);
1678   return res;
1679 }
1680
1681 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)
1682 {
1683 // create a 2D histogram
1684
1685   TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
1686   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1687   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1688   //  res->SetMarkerStyle(kFullCircle);
1689   //  res->SetOption("E");
1690   res->SetLineColor(color);
1691   //  fOutputList->Add(res);
1692   return res;
1693 }
1694
1695
1696 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)
1697 {
1698 // create a 2D histogram
1699
1700   TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
1701   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1702   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1703   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1704   //  res->SetMarkerStyle(kFullCircle);
1705   //  res->SetOption("E");
1706   res->SetLineColor(color);
1707   //  fOutputList->Add(res);
1708   return res;
1709 }
1710
1711 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)
1712 {
1713 // create a 2D histogram
1714
1715   TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
1716   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1717   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1718   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1719   //  res->SetMarkerStyle(kFullCircle);
1720   //  res->SetOption("E");
1721   res->SetLineColor(color);
1722   //  fOutputList->Add(res);
1723   return res;
1724 }
1725
1726 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1727                                                   Int_t nBinsy, Double_t yMin, Double_t yMax,
1728                                                   Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1729 {
1730 // create a 3D histogram
1731
1732   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);
1733   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1734   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1735   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1736   //  res->SetMarkerStyle(kFullCircle);
1737   //  res->SetOption("E");
1738   res->SetLineColor(color);
1739   //  fOutputList->Add(res);
1740   return res;
1741 }
1742
1743 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1744                                                   Int_t nBinsy, Double_t yMin, Double_t yMax,
1745                                                   Int_t nBinsz, Double_t *zArr, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1746 {
1747 // create a 3D histogram
1748
1749   Double_t xArr[nBinsx+1], yArr[nBinsy+1];
1750
1751   for(Int_t i=0; i<=nBinsx; i++) xArr[i]=xMin+i*(xMax-xMin)/nBinsx;
1752   for(Int_t i=0; i<=nBinsy; i++) yArr[i]=yMin+i*(yMax-yMin)/nBinsy;
1753
1754   TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xArr, nBinsy, yArr, nBinsz, zArr);
1755   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1756   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1757   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1758   //  res->SetMarkerStyle(kFullCircle);
1759   //  res->SetOption("E");
1760   res->SetLineColor(color);
1761   //  fOutputList->Add(res);
1762   return res;
1763 }
1764
1765
1766 //////////////////////////////////////
1767
1768
1769 /*
1770 //________________________________________________________________________
1771 void AnalysisJetMain::ConnectInputData(Option_t *)
1772 {
1773   // Connect ESD
1774   // Called once
1775
1776   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1777   //   if (!fESD) {
1778   //     AliError("ESD not available");
1779   //     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1780   //     fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1781
1782   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1783   //  assume that the AOD is in the general output...
1784   fAODOut  = AODEvent();
1785
1786
1787 }
1788 */
1789
1790 void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1791 {
1792    fJetBranchName[0] = branch1;
1793    fJetBranchName[1] = branch2;
1794 }
1795
1796 //________________________________________________________________________
1797 void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1798 {
1799   //printf("Open1\n");
1800   //       const char *nameF = OpenFile(1)->GetName();
1801   //printf("Open2 %s\n",nameF);
1802
1803   //  fTriggerAnalysis = new AliTriggerAnalysis();
1804   //  if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1805
1806
1807   fOutputList = new TList();
1808   fOutputList->SetOwner();
1809
1810   fhPtJL  = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1811   fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1812
1813
1814
1815   printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1816
1817
1818   fanalJetSubStrHM->SetFilterMask(fFilterMask);
1819   if(fkMC)   fanalJetSubStrHM->MCprod();
1820   fanalJetSubStrHM->InitHistos();
1821   fanalJetSubStrHM->AddToList(fOutputList);
1822
1823
1824   if(fkMC)
1825     {
1826
1827       printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1828       fanalJetSubStrMCHM->InitHistos();
1829       fanalJetSubStrMCHM->AddToList(fOutputList);
1830
1831       for(Int_t i=0; i<3; i++)
1832         {
1833           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]);
1834           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]);
1835           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]);
1836           fhDCAxy[i]        = Hist2D(Form("hDCAxy%d",i), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]"); fOutputList->Add(fhDCAxy[i]);
1837           fhDCAz[i]         = Hist2D(Form("hDCAz%d",i) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ; fOutputList->Add(fhDCAz[i]);
1838           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]);
1839         }
1840     }
1841
1842
1843
1844
1845    printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1846
1847
1848
1849   PostData(1, fOutputList);
1850
1851   /*
1852    OpenFile(1);
1853
1854   fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1855   fOutputTree->Branch("event",&f2JEvent);
1856   */
1857
1858 //    PostData(1, fOutputTree);
1859 }
1860
1861 /*
1862 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1863 {
1864
1865   // Fetch the aod also from the input in,
1866   // have todo it in notify
1867
1868
1869   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1870   //  assume that the AOD is in the general output...
1871   fAODOut  = AODEvent();
1872
1873   if(fNonStdFile.Length()!=0){
1874     // case that we have an AOD extension we need can fetch the jets from the extended output
1875     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1876     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1877     if(!fAODExtension){
1878       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1879     }
1880   }
1881
1882
1883
1884   return kTRUE;
1885 }
1886 */
1887
1888
1889
1890
1891 //________________________________________________________________________
1892 void AliAnalysisTaskJetShape::UserExec(Option_t *)
1893 {
1894   //   return;
1895 //    if(f2JEvent)
1896 //      delete f2JEvent;
1897 //    f2JEvent = new Ali2JEvent();
1898 //    return;
1899
1900   //  return;
1901  if(fDebug)  Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
1902   if(fDebug) printf("NEW EVENT 0\n");
1903
1904   if(!IsGoodEvent()) return;
1905
1906
1907   //   printf("\n\n\n NEW EVENT\n");
1908
1909    AliAODEvent* aodE = 0;
1910    if(!aodE){
1911      if(!fESD)aodE = fAODIn;
1912      else aodE = fAODOut;}
1913
1914    /*
1915    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1916 if(!aodH){
1917     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1918     return;
1919   }
1920    */
1921    if(fDebug) {
1922      printf("NEW EVENT 1:  Number of tracks = %d\n",aodE->GetNumberOfTracks());
1923      //     aodE->GetList()->Print();
1924      //     printf("Look for %s\n",fJetBranchName[0].Data());
1925    }
1926
1927
1928   // centrality selection
1929    AliCentrality *cent = 0x0;
1930    Double_t centrality = 0.;
1931
1932    if(fESD) {cent = fESD->GetCentrality();
1933      if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1934    else     centrality=((AliVAODHeader*)aodE->GetHeader())->GetCentrality();
1935
1936
1937    if(!fkIsPbPb) {
1938      centrality=1.;
1939    }
1940
1941 //  if(fDebug) printf("NEW EVENT 2\n");
1942
1943    Bool_t  fUseAOD049 = kFALSE;// kTRUE;//
1944       if(fUseAOD049&&centrality>=0){
1945         Float_t v0=0;
1946         AliAODVZERO* aodV0 = aodE->GetVZEROData();
1947         v0+=aodV0->GetMTotV0A();
1948         v0+=aodV0->GetMTotV0C();
1949         if(centrality==0 && v0 < 19500) return ;//filtering issue
1950         Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
1951         Float_t val= 1.30552 +  0.147931 * v0;
1952         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,
1953             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,
1954             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,
1955             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,
1956             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,
1957             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,
1958             13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1959         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] )
1960           return; //outlier
1961       }
1962
1963
1964     if(fDebug)
1965       {
1966         printf("centrality: %f\n", centrality);
1967         aodE->Print();
1968         aodE->GetList()->Print();
1969       }
1970
1971
1972     if (centrality < fCentMin || centrality > fCentMax){
1973     //    PostData(1, fOutputList);
1974     return;
1975     }
1976
1977    //   fhNEvent->Fill(0);
1978
1979    // accepted events
1980    // -- end event selection --
1981
1982
1983    // get background
1984    AliAODJetEventBackground* externalBackground = 0;
1985    AliAODJetEventBackground* externalBackgroundMC = 0;
1986
1987
1988    Float_t rho = 0;
1989    Float_t rho_MC=0;
1990
1991    if(fkIsPbPb)
1992      {
1993        if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
1994          externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
1995          if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1996        }
1997        if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
1998          externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
1999          if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2000        }
2001        if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
2002          externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
2003          if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2004        }
2005        if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2006          externalBackgroundMC =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
2007          if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2008        }
2009        if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2010          externalBackgroundMC =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
2011          if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2012        }
2013        if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2014          externalBackgroundMC =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
2015          if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2016        }
2017        if(externalBackground)rho = externalBackground->GetBackground(0);
2018        if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
2019      }
2020
2021
2022
2023   if(fkMC)
2024     {
2025       AliAODVertex* primVtx = aodE->GetPrimaryVertex();
2026       Double_t bfield = aodE->GetMagneticField();
2027
2028       TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2029       if(!arrP)
2030         {
2031           printf("ERROR: no Info about particles in AOD\n");
2032           return;
2033         }
2034
2035       for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2036         {
2037           AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodE->GetTrack(it));
2038           if(!tr) AliFatal("Not a standard AOD");
2039           if(!tr) continue;
2040           if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2041           if(TMath::Abs(tr->Eta())>1.) continue;
2042           Int_t label = TMath::Abs(tr->GetLabel());
2043
2044           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2045           if(!part)continue;
2046           //      if(!part->IsPhysicalPrimary())continue;
2047           //      if(part->Charge()==0)continue;
2048           //      vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2049
2050
2051
2052           Double_t dca[2];
2053           Double_t cov[3];
2054
2055           AliAODTrack tmp(*tr);
2056           tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2057
2058           Int_t type = 0;
2059           if(!part->IsPhysicalPrimary()) type=1;
2060           if(label<0) type =2;
2061
2062           fhPtresVsPt[type]->Fill(part->Pt(), 1-tr->Pt()/part->Pt());
2063           fhPhiresVsPhi[type]->Fill(part->Phi(), tr->Phi() - part->Phi());
2064           fhEtaresVsEta[type]->Fill(part->Eta(), tr->Eta() - part->Eta());
2065           fhDCAxy[type]->Fill(part->Pt(), dca[0]);
2066           fhDCAz[type]->Fill( part->Pt(), dca[1]);
2067           fhTrackPtEtaPhi[type]->Fill(tr->Pt(), tr->Eta(), tr->Phi());
2068         }
2069
2070
2071     }
2072
2073
2074
2075
2076
2077    //   printf("rho = %f\n",rho);
2078
2079     //   return;
2080
2081    TClonesArray *Jets[2];
2082    Jets[0]=0;
2083    Jets[1]=0;
2084    if(fAODOut&&!Jets[0]){
2085    Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
2086    Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
2087    if(fAODExtension && !Jets[0]){
2088    Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
2089    Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
2090    if(fAODIn&&!Jets[0]){
2091    Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
2092    Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
2093
2094
2095    if(!Jets[0]) {
2096      Printf("no friend!!!\n");
2097      return;
2098    }
2099   Int_t nJ = Jets[0]->GetEntries();
2100   Float_t ptmax = 0.;
2101   Int_t   imax  = -1;
2102
2103   if(fDebug) {
2104     printf("NEW EVENT 3:  Number of Rec. jets  %d\n",nJ);
2105   }
2106
2107 //
2108 // Find highest pT jet with pt > 20 GeV
2109 //
2110 //  fPtJcorrMin=0;
2111
2112   Float_t etaJmax = 0.4;
2113
2114      for (Int_t i = 0; i < nJ; i++) {
2115        AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2116        if (!jet) continue;
2117        //              jet->Print("0");
2118        Float_t ptJ  = jet->Pt();
2119        Float_t etaJ = TMath::Abs(jet->Eta());
2120
2121        Float_t area = jet->EffectiveAreaCharged();
2122        Float_t ptJcorr=ptJ-rho*area;
2123
2124        if ((ptJcorr > fPtJcorrMin) && (ptJcorr  > ptmax) && etaJ < etaJmax) {
2125          ptmax = ptJcorr;
2126          imax = i;
2127        }
2128
2129      }
2130
2131
2132   if(fDebug) {
2133     printf("NEW EVENT 4:\n");
2134   }
2135
2136      TVector3 vecJ;
2137
2138      AliAODJet* jetL = 0;
2139
2140      Float_t areaJL, ptJLcorr;
2141
2142      if (imax != -1)  {
2143
2144        jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2145        if(jetL){
2146
2147          areaJL   = jetL->EffectiveAreaCharged();
2148          ptJLcorr  = jetL->Pt()-rho*areaJL;
2149
2150          fhPtJL->Fill(ptJLcorr);
2151          fhAreaJL->Fill(areaJL);
2152          vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2153          fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
2154
2155          if(fDebug) {
2156            Printf("Leading Jet");
2157            jetL->Print("0");
2158          }
2159
2160        }
2161      }
2162
2163
2164      if(!fkMC)
2165        {
2166          PostData(1, fOutputList);
2167          if(fDebug)  Printf("End of event #%5d", (Int_t) fEntry);
2168          return;
2169        }
2170
2171
2172
2173
2174      nJ = Jets[1]->GetEntries();
2175   if(fDebug) {
2176     printf("NEW EVENT 5:  Number of Rec. jets  %d\n",nJ);
2177   }
2178
2179      ptmax = 0;
2180      imax = -1;
2181      Float_t areaJL_MC=0;
2182
2183      for (Int_t i = 0; i < nJ; i++) {
2184        AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2185        if (!jet) continue;
2186        Float_t ptJ1  = jet->Pt();
2187        Float_t etaJ1 = TMath::Abs(jet->Eta());
2188
2189        areaJL_MC = jet->EffectiveAreaCharged();
2190        Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2191
2192        //    jet->Print();
2193
2194        if ((ptJcorr > fPtJcorrMin) && (ptJcorr  > ptmax) && etaJ1 < etaJmax) {
2195          ptmax = ptJcorr;
2196          imax = i;
2197        }
2198
2199      }
2200
2201   if(fDebug) {
2202     printf("NEW EVENT 6:\n");
2203   }
2204
2205
2206      AliAODJet* jetMCL=0;
2207
2208     if (imax != -1)  {
2209       jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
2210       if(jetMCL){
2211         fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2212       }
2213     }
2214
2215
2216
2217
2218
2219      PostData(1, fOutputList);
2220
2221     if(fDebug)  Printf("End of event #%5d", (Int_t) fEntry);
2222
2223      return;
2224
2225
2226
2227
2228 }
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242 Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2243 {
2244   while(phi1<0) phi1+=TMath::TwoPi();
2245   while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2246
2247   while(phi2<0) phi2+=TMath::TwoPi();
2248   while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
2249
2250   Double_t dphi = phi1- phi2;
2251   if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
2252   if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
2253
2254   //  Double_t dphi = phi1- phi2;
2255   //  while(dphi<0) dphi+=TMath::TwoPi();
2256   //  while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2257   return dphi;
2258 }
2259
2260 Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2261 {
2262
2263   Double_t dphi = CalcdPhi(phi1, phi2);
2264   while(dphi<0) dphi+=TMath::TwoPi();
2265   while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2266   return dphi;
2267 }
2268
2269
2270 Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2271 {
2272
2273   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2274   //   if (!fESD) {
2275   //     AliError("ESD not available");
2276   //     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2277   //     fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2278
2279   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2280   //  assume that the AOD is in the general output...
2281   fAODOut  = AODEvent();
2282
2283
2284    static AliAODEvent* aod = 0;
2285    // take all other information from the aod we take the tracks from
2286    if(!aod){
2287      if(!fESD)aod = fAODIn;
2288      else aod = fAODOut;}
2289
2290
2291
2292    if(fNonStdFile.Length()!=0){
2293     // case that we have an AOD extension we need can fetch the jets from the extended output
2294      AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2295      fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2296    }
2297
2298    if(fDebug){
2299     if(fAODIn) printf("fAODIn\n");
2300     if(fAODOut) printf("fAODOut\n");
2301     if(fAODExtension) printf("fAODExtension\n");
2302    }
2303
2304
2305    /*
2306    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2307 if(!aodH){
2308     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2309       return kFALSE;
2310   }
2311    */
2312
2313    // -- event selection --
2314
2315    // physics selection
2316    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2317    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2318    //     cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
2319    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2320       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2321       return kFALSE;
2322    }
2323
2324    // vertex selection
2325    //   if(!aod){
2326    //     if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2327    //     fhNEvent->Fill(3);
2328    //      PostData(1, fOutputList);
2329    //     return kFALSE;
2330    //   }
2331
2332    AliAODVertex* primVtx = aod->GetPrimaryVertex();
2333
2334    if(!primVtx){
2335      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
2336      return kFALSE;
2337    }
2338
2339    Int_t nTracksPrim = primVtx->GetNContributors();
2340    if ((nTracksPrim < fVtxMinContrib) ||
2341          (primVtx->GetZ() < fVtxZMin) ||
2342          (primVtx->GetZ() > fVtxZMax) ){
2343       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2344       return  kFALSE;
2345    }
2346
2347    /*
2348    // event class selection (from jet helper task)
2349    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
2350    if(fDebug) Printf("Event class %d", eventClass);
2351    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
2352       return  kFALSE;
2353    }
2354    */
2355     return kTRUE;
2356 }
2357
2358
2359
2360
2361
2362
2363 //________________________________________________________________________
2364 void AliAnalysisTaskJetShape::Terminate(Option_t *)
2365 {
2366   // Draw result to the screen
2367   // Called once at the end of the query
2368  //    fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2369
2370   printf("MyTaskTestTrigger::Terminate()\n\n\n");
2371
2372
2373   //  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2374   //  fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2375
2376 }
2377
2378
2379
2380 /*
2381 Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd,  Double_t Vxyz[3], Int_t type)
2382 {
2383   // type
2384   //     0 -> SPD vtx
2385   //     1 -> ESD vtx
2386   const AliESDVertex* vtx = 0;
2387
2388   if(type==0) {
2389     vtx = esd->GetPrimaryVertexSPD();
2390     if(!vtx) return kFALSE;
2391     if(vtx->GetNContributors() < 1) return kFALSE;
2392     TString vtxTyp = vtx->GetTitle();
2393     if (vtxTyp.Contains("vertexer: Z")) {
2394       if (vtx->GetDispersion()>0.04) return kFALSE;
2395       if (vtx->GetZRes()>0.25) return kFALSE;
2396     }
2397
2398   }
2399   if(type==1) {
2400     vtx = esd->GetPrimaryVertexTracks();
2401     if(!vtx) return kFALSE;
2402     if(vtx->GetNContributors() < 1) return kFALSE;
2403   }
2404
2405
2406    Vxyz[0] = vtx->GetX();
2407    Vxyz[1] = vtx->GetY();
2408    Vxyz[2] = vtx->GetZ();
2409
2410    return kTRUE;
2411 }
2412 */
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422 TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax,  const char* xLabel, Int_t color)
2423 {
2424 // create a 1D histogram
2425
2426   TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
2427   if (xLabel) h->GetXaxis()->SetTitle(xLabel);
2428   //  res->SetMarkerStyle(kFullCircle);
2429   //  res->SetOption("E");
2430   h->SetLineColor(color);
2431   //  fOutputList->Add(h);
2432   return h;
2433 }
2434
2435
2436 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)
2437 {
2438 // create a 2D histogram
2439
2440   TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
2441   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2442   if (xLabel) res->GetYaxis()->SetTitle(yLabel);
2443   //  res->SetMarkerStyle(kFullCircle);
2444   //  res->SetOption("E");
2445   res->SetLineColor(color);
2446   //  fOutputList->Add(res);
2447   return res;
2448 }
2449
2450
2451 TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
2452                                                   Int_t nBinsy, Double_t yMin, Double_t yMax,
2453                                                   Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
2454 {
2455 // create a 3D histogram
2456
2457   TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
2458   if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2459   if (yLabel) res->GetYaxis()->SetTitle(yLabel);
2460   if (zLabel) res->GetZaxis()->SetTitle(zLabel);
2461   //  res->SetMarkerStyle(kFullCircle);
2462   //  res->SetOption("E");
2463   res->SetLineColor(color);
2464   //  fOutputList->Add(res);
2465   return res;
2466 }
2467
2468
2469
2470
2471 //#endif
2472
2473
2474
2475
2476
2477
2478
2479