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