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