Adding Jet Shape User task (M. Poghosyan)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Nov 2012 13:29:49 +0000 (13:29 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Nov 2012 13:29:49 +0000 (13:29 +0000)
PWGJE/CMakelibPWGJE.pkg
PWGJE/PWGJELinkDef.h
PWGJE/UserTasks/AliAnalysisTaskJetShape.cxx [new file with mode: 0644]
PWGJE/UserTasks/AliAnalysisTaskJetShape.h [new file with mode: 0644]

index 98235c5..946e0e8 100644 (file)
@@ -47,6 +47,7 @@ set ( SRCS
     UserTasks/AliAnalysisTaskCheckSingleTrackJetRejection.cxx 
     UserTasks/AliAnalysisTaskJetHadronCorrelation.cxx
     UserTasks/AliAnalysisTaskJetHBOM.cxx
+    UserTasks/AliAnalysisTaskJetShape.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
index c090e86..161cce9 100644 (file)
@@ -33,4 +33,7 @@
 #pragma link C++ class AliAnalysisTaskCheckSingleTrackJetRejection+;
 #pragma link C++ class AliAnalysisTaskJetHadronCorrelation+;
 #pragma link C++ class AliAnalysisTaskJetHBOM+;
+#pragma link C++ class AliAnalysisTaskJetShape+;
+#pragma link C++ class AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool+;
+#pragma link C++ class AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM+;
 #endif
diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetShape.cxx b/PWGJE/UserTasks/AliAnalysisTaskJetShape.cxx
new file mode 100644 (file)
index 0000000..5a50981
--- /dev/null
@@ -0,0 +1,2134 @@
+// Class AliAnalysisTaskJetShape
+// Author martin.poghosyan@cern.ch
+//
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TCanvas.h"
+#include "TParticlePDG.h"
+#include "TParticle.h"
+#include "TRandom.h"
+#include "TDatabasePDG.h"
+#include <TPDGCode.h>
+#include "THnSparse.h"
+#include "TGraph.h" 
+#include "TArrayF.h"
+#include "TArrayD.h"
+#include <TVector3.h>
+#include <TClonesArray.h>
+#include "AliLog.h"
+
+#include "AliCDBManager.h"
+#include "AliGeomManager.h"
+#include "AliCDBEntry.h"
+#include "AliCDBPath.h"
+#include "AliAlignObjParams.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliMultiplicity.h"
+#include "AliESDFMD.h"
+#include "AliESDVZERO.h"
+#include "AliESDInputHandler.h"
+#include "AliVEvent.h"
+#include "AliCentrality.h"
+#include "AliAnalysisHelperJetTasks.h"
+#include "AliInputEventHandler.h"
+#include "AliAODJetEventBackground.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODJet.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliGenEventHeader.h"
+#include "AliHeader.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+#include "AliInputEventHandler.h"
+#include "AliLHCData.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliESDtrackCuts.h"
+#include "AliTriggerAnalysis.h"
+
+
+
+#include "AliAnalysisTaskJetShape.h"
+
+/////////////////////////////////////
+
+ClassImp(AliAnalysisTaskJetShape)
+
+
+  AliAnalysisTaskJetShape::AliAnalysisTaskJetShape(const char *name) 
+: AliAnalysisTaskSE(name), 
+  //  fOutputTree(0),
+  fOutputList(0),
+  fESD(0x0),
+  fAODIn(0x0),
+  fAODOut(0x0),
+  fAODExtension(0x0),
+  fkMC(0),       
+  fCMSE(0),      
+  fRunNb(0),       
+  fkIsPhysSel(0),  
+  fNonStdFile(""),
+  fFilterMask(0),
+  fOfflineTrgMask(AliVEvent::kAny),
+  fCentMin(0.),
+  fCentMax(80.),
+  fEvtClassMin(0),
+  fEvtClassMax(4),
+  fPtJcorrMin(20),
+  fPtJBcorrMin(20),
+  fJpPtmin(1),
+  fJaPtmin(0.5),
+  fVtxMinContrib(1),
+  fVtxZMin(-10),
+  fVtxZMax(10),
+  fkIsPbPb(0),
+  fhPtJL(0),
+  fhAreaJL(0),
+  fanalJetSubStrHM(0),
+  fanalJetSubStrMCHM(0)
+{
+
+
+  fgkbinCent[0] = -0.099;
+  fgkbinCent[1] = 5;
+  fgkbinCent[2] = 10;
+  fgkbinCent[3] = 20;
+  fgkbinCent[4] = 40;
+  fgkbinCent[5] = 60;
+  fgkbinCent[6] = 80;
+  fgkbinCent[7] =100;
+
+  fBackgroundBranch[0] = "",
+  fBackgroundBranch[1] = "",
+  fJetBranchName[0] = "";
+  fJetBranchName[1] = "";
+
+  //  fkIsPbPb = kFALSE;
+  //  fDebug=0;
+
+
+
+    fanalJetSubStrHM   = new AliAnalysisTaskJetShapeHM("rec");
+    fanalJetSubStrMCHM = new AliAnalysisTaskJetShapeHM("truth", kTRUE);
+
+  DefineOutput(1, TList::Class());
+
+  //   DefineOutput(1, TTree::Class());
+}
+
+
+AliAnalysisTaskJetShape::~AliAnalysisTaskJetShape()
+{
+   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+    printf("Deleteing output\n");
+
+    if(fOutputList){
+      delete fOutputList;
+      fOutputList = 0;
+    }
+
+    //    if(fTriggerAnalysis) 
+    //      delete fTriggerAnalysis;
+
+   } 
+}
+
+
+AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool():TObject(),
+  fvecJ(0,0,0),
+  fvecB1(0,0,0),
+  fvecB2(0,0,0),
+  fRmax(0),
+  fPRecInRJ(0,0,0),
+  fList(0),
+  fEntries(0)
+{
+  fList = new TClonesArray("TVector3",10000);
+
+  TVector3 v(0,0,0);
+  SetVecJ(v);
+  fRmax = -0.5;
+  fR[0] =0.1;
+  fR[1] =0.2;
+  fR[2] =0.3;
+  fEntries=0;
+
+  for(Int_t i1=0; i1<fgkbinR; i1++) 
+    {
+      for(Int_t i2=0; i2<2; i2++)
+       {
+         fListJ[i1][i2].Set(1000); 
+         fListB1[i1][i2].Set(1000);
+         fListB2[i1][i2].Set(1000);
+         fListJc[i1][i2]  = 0; 
+         fListB1c[i1][i2] = 0;
+         fListB2c[i1][i2] = 0;
+         //      fListJ[i1][i2]("TVector3",10000); 
+         //      fListB1[i1][i2] = TClonesArray("TVector3",10000);
+         //      fListB2[i1][i2] = TClonesArray("TVector3",10000);
+
+                 fListJ[i1][i2].Reset(-1); 
+                 fListB1[i1][i2].Reset(-1);
+                 fListB2[i1][i2].Reset(-1);
+
+                 fPRecJ[i1][i2].SetXYZ(0,0,0);
+
+       }
+    }
+
+  fPRecInRJ.SetXYZ(0,0,0);
+
+  fPtMinTr[0] = 2;
+  fPtMinTr[1] = 0.5;
+
+
+}
+
+AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool(TClonesArray *list):TObject(),
+  fvecJ(0,0,0),
+  fvecB1(0,0,0),
+  fvecB2(0,0,0),
+  fRmax(0),
+  fPRecInRJ(0,0,0),
+  fList(list),
+  fEntries(0)
+{
+
+  TVector3 v(0,0,0);
+  SetVecJ(v);
+  fRmax = 0.4;
+  fR[0] =0.1;
+  fR[1] =0.2;
+  fR[2] =0.3;
+  fEntries=0;
+
+  for(Int_t i1=0; i1<fgkbinR; i1++) 
+    {
+      for(Int_t i2=0; i2<2; i2++)
+       {
+         fListJ[i1][i2].Set(1000); 
+         fListB1[i1][i2].Set(1000);
+         fListB2[i1][i2].Set(1000);
+         fListJc[i1][i2]  = 0; 
+         fListB1c[i1][i2] = 0;
+         fListB2c[i1][i2] = 0;
+         //      fListJ[i1][i2]("TVector3",10000); 
+         //      fListB1[i1][i2] = TClonesArray("TVector3",10000);
+         //      fListB2[i1][i2] = TClonesArray("TVector3",10000);
+                 fListJ[i1][i2].Reset(-1); 
+                 fListB1[i1][i2].Reset(-1);
+                 fListB2[i1][i2].Reset(-1);
+                 fPRecJ[i1][i2].SetXYZ(0,0,0);
+
+       }
+    }
+
+  fPRecInRJ.SetXYZ(0,0,0);
+
+  fPtMinTr[0] = 2;
+  fPtMinTr[1] = 0.5;
+
+}
+
+
+
+
+AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::~AliAnalysisTaskJetShapeTool()
+{
+  fList->Delete();
+  fEntries=0;
+
+  //  if(fList)
+  //    delete fList;
+
+}
+
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::SetVecJ(TVector3 v)
+{
+  fvecJ.SetXYZ(v.X(), v.Y(), v.Z());
+  fvecB1.SetXYZ(-v.Y(), v.X(), v.Z());
+  fvecB2.SetXYZ(v.Y(),-v.X(), v.Z());
+}
+
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Clean()
+{
+  //  printf("AnalJetSubStrTool::Clean()\n");
+  /*
+  fList->Delete();
+  for(Int_t i1=0; i1<fgkbinR; i1++) 
+    {
+      for(Int_t i2=0; i2<2; i2++)
+       {
+         fListJ[i1][i2]->Delete(); 
+         fListB1[i1][i2]->Delete();
+         fListB2[i1][i2]->Delete();
+       }
+    }
+
+  */
+
+  //  fList.SetOwner();
+  //  fList->Delete();
+  //  fEntries=0;
+  fPRecInRJ.SetXYZ(0,0,0);
+
+  for(Int_t i1=0; i1<fgkbinR; i1++) 
+    {
+      for(Int_t i2=0; i2<2; i2++)
+       {
+                 fListJ[i1][i2].Reset(-1); 
+                 fListB1[i1][i2].Reset(-1);
+                 fListB2[i1][i2].Reset(-1);
+         //      fListJ[i1][i2].Set(0); 
+         //      fListB1[i1][i2].Set(0);
+         //      fListB2[i1][i2].Set(0);
+         fListJc[i1][i2]  = 0; 
+         fListB1c[i1][i2] = 0;
+         fListB2c[i1][i2] = 0;
+
+         fPRecJ[i1][i2].SetXYZ(0,0,0);
+       }
+    }
+  
+}
+
+
+//void AnalJetSubStrTool::Print(Option_t* /*option*/) const
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT() const
+{
+
+  //     PRINT(fList, "all"); 
+
+         //  fList->Print();
+  for(Int_t i1=0; i1<fgkbinR; i1++) 
+    {
+      if(i1!=0) continue;
+     for(Int_t i2=0; i2<2; i2++)
+       {
+
+         //      printf("# %d   %d\n",i1, i2);
+         PRINT(fListJ[i1][i2], fListJc[i1][i2], Form("J_%d_%d",i1,i2)); 
+                 //              PRINT(fListB1[i1][i2], Form("B1_%d_%d",i1,i2)); 
+                 //              PRINT(fListB2[i1][i2], Form("B2_%d_%d",i1,i2)); 
+         //      fListJ[i1][i2].Print("",1); 
+         //      fListB1[i1][i2]->Print();
+         //      fListB2[i1][i2]->Print();
+       }
+    }
+
+}
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT(TArrayI a, Int_t n, const char *txt) const
+{
+  printf("%s :%d \n",txt, n);
+  Int_t count = 0;
+  for(Int_t i=0; i<n; i++){
+    printf("%4d   ", a.At(i));
+
+    if(count==20) 
+      {
+       printf("\n");
+       count=0;
+      }
+    else count++;
+  }
+  if(count!=20) printf("\n");
+
+}
+
+
+
+ void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Init()
+{
+  Int_t Nev = fEntries;
+
+  //    PRINT();
+
+  //  printf("Nev = %d\n",Nev);
+
+  for(Int_t iP=0; iP<Nev; iP++) 
+    {
+      TVector3 *v = (TVector3*) fList->At(iP);
+      if(!v) {
+       printf("ERROR : entry #%d not found\n",iP);
+       continue;
+      }
+
+      //      printf("#%3d   ",iP); v->Print();
+
+      Double_t R = CalcR(fvecJ, *v);
+      //     printf("R = %f\n",R);
+      if(R<fRmax)
+       {
+         for(Int_t iT = 0; iT < fgkbinR; iT++)
+           {
+             Int_t type = 0;
+             if(R>fR[iT]) type = 1;
+
+             if(v->Pt() < fPtMinTr[type]) continue;
+             fPRecJ[iT][type]+=*v;
+             AddToJ(iT, type, iP);
+           }
+
+         if(v->Pt() < fPtMinTr[0]) continue;
+         fPRecInRJ+=*v;
+
+         continue;
+       }
+
+      R = CalcR(fvecB1, *v);
+      if(R<fRmax)
+       {
+         for(Int_t iT = 0; iT < fgkbinR; iT++)
+           {
+             Int_t type = 0;
+             if(R>fR[iT]) type = 1;
+
+             if(v->Pt() < fPtMinTr[type]) continue;
+
+             AddToB1(iT, type, iP);
+           }
+         continue;
+       }
+
+      R = CalcR(fvecB2, *v);
+      if(R<fRmax)
+       {
+         for(Int_t iT = 0; iT < fgkbinR; iT++)
+           {
+             Int_t type = 0;
+             if(R>fR[iT]) type = 1;
+
+             if(v->Pt() < fPtMinTr[type]) continue;
+
+             AddToB2(iT, type, iP);
+           }
+         continue;
+       }
+    }
+
+
+
+  /*
+  for(Int_t i1=0; i1<fgkbinR; i1++) 
+    {
+      for(Int_t i2=0; i2<2; i2++)
+       {
+         //      if(fListJc[i1][i2]<2) fListJc[i1][i2]=0;
+         //      if(fListB1c[i1][i2]<2) fListB1c[i1][i2]=0;
+         //      if(fListB2c[i1][i2]<2) fListB2c[i1][i2]=0;
+         fListJ[i1][i2].Set(fListJc[i1][i2]); 
+         fListB1[i1][i2].Set(fListB1c[i1][i2]);
+         fListB2[i1][i2].Set(fListB2c[i1][i2]);
+       }
+    }
+  */
+
+}
+
+
+
+Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcR(TVector3 v1, TVector3 v2)
+{
+
+  Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
+  //  dphi*=(0.9/TMath::Pi());
+  Double_t deta = v1.Eta() - v2.Eta();
+  Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
+  return RB;
+}
+
+
+Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcdPhi(Double_t phi1, Double_t phi2)
+{
+
+  while(phi1<0) phi1+=TMath::TwoPi();
+  while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
+
+  while(phi2<0) phi2+=TMath::TwoPi();
+  while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
+
+  Double_t dphi = phi1- phi2;
+  if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
+  if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
+
+  return dphi;
+}
+
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::GetLocalMom(TVector3 vecJ, TVector3 *q)
+{
+  // theta and phi of a particle in loc.syst of fvecJ 
+
+
+  Double_t PT = vecJ.Perp();
+  Double_t P0 = vecJ.Mag();
+
+  Double_t q0=q->X();
+  Double_t q1=q->Y();
+  Double_t q2=q->Z();
+
+  Double_t qx1 = vecJ(2)/P0/PT*(vecJ(0)*q0+vecJ(1)*q1) - PT/P0*q2;
+  Double_t qy1 = -vecJ(1)/PT*q0 + vecJ(0)/PT*q1;
+  Double_t qz1 = 1/P0*(vecJ(0)*q0+vecJ(1)*q1+vecJ(2)*q2);
+
+  //  Double_t q00 = TMath::Sqrt(qx1*qx1 + qy1*qy1 + qz1*qz1);
+  //  phi = TMath::ATan2(qy1, qx1);
+  //  if(phi<0) phi+=fTwoPi;
+  //  theta = TMath::ACos(qz1/q00);
+
+  q->SetXYZ(qx1, qy1, qz1);
+
+  return;
+
+}
+
+/*
+
+Bool_t AnalJetSubStrTool::FindCorrelationAxesAnd(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
+{
+  Double_t TwoPi = TMath::TwoPi();
+
+//
+// 1st track loop to determine the sphericity matrix
+//   
+      // Initialize Shericity matrix
+      Float_t mxx    = 0.;
+      Float_t myy    = 0.;
+      Float_t mxy    = 0.;
+      Int_t   nc     = 0;
+      Float_t sump2  = 0.;
+      Float_t ptMax  = 0.;
+      Float_t etaMax = 0.;
+      Float_t phiMax = 0.;
+      Int_t   iMax   = -1;
+      Float_t ptsum  = 0.;
+
+  Int_t Nev =  list.GetSize();
+  if( Nev <2) return kFALSE;
+
+  Int_t indexpmax = -1;
+  Double_t pmax = -1;
+  Double_t phipmax = 0;
+
+  Int_t indexpzmax = -1;
+  Double_t pzmax = -1;
+  Double_t phipzmax = 0;
+
+
+  Int_t indexthetamin = -1;
+  Double_t thetamin = 1000;
+  Double_t phithetamin = 0;
+
+  for(Int_t ip=0; ip< Nev; ip++) {
+
+    TVector3* part = (TVector3*) GetAt(list.At(ip));
+    if(!part) continue;
+
+    TVector3 vtmp(*part);
+    GetLocalMom(vec, &vtmp);
+
+    Float_t ppjX = vtmp.X();
+    Float_t ppjY = vtmp.Y();
+    Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
+
+
+    Float_t pt = ppjT;//part->Z();
+    Float_t deta = part->Eta();
+    Float_t dphi = part->Phi();
+
+             mxx += (ppjX * ppjX / ppjT);
+             myy += (ppjY * ppjY / ppjT);
+             mxy += (ppjX * ppjY / ppjT);
+             nc++;
+             sump2 += ppjT;
+
+
+    if(vtmp.Mag()>pmax)
+      {
+       pmax=vtmp.Mag();
+       indexpmax = ip;
+       phipmax=vtmp.Phi();
+      }
+
+    if(vtmp.Z()>pzmax)
+      {
+       pzmax=vtmp.Z();
+       indexpzmax = ip;
+       phipzmax=vtmp.Phi();
+      }
+    if(vtmp.Theta()<thetamin)
+      {
+       thetamin=vtmp.Theta();
+       indexthetamin = ip;
+       phithetamin=vtmp.Phi();
+      }
+
+
+         }
+  
+
+
+    
+//
+// At this point we have mxx, myy, mxy
+  if (nc == 0) return kFALSE;      
+// Shericity Matrix    
+      const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};    
+      TMatrixDSym m0(2, ele);
+// Find eigenvectors
+      TMatrixDSymEigen m(m0);
+      TVectorD eval(2);
+      TMatrixD evecm = m.GetEigenVectors();
+      eval  = m.GetEigenValues();
+// Largest eigenvector
+      Int_t jev = 0;
+      if (eval[0] < eval[1]) jev = 1;
+      TVectorD evec0(2);
+// Principle axis
+      evec0 = TMatrixDColumn(evecm, jev);
+      TVector2 evec(evec0[0], evec0[1]); 
+// Principle axis from leading partice
+
+//      Float_t phiM = TMath::ATan2(phiMax, etaMax);
+//      TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM)); 
+//      Float_t phistM = evecM.DeltaPhi(evec);
+//      if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
+
+
+      Double_t phiTA = evec.Phi();
+
+
+      Double_t phiDir = 0;
+  if(scenario ==0)
+    {
+      phiDir = phipmax;
+    }
+  else if(scenario ==1)
+    {
+      phiDir = phipzmax;
+    }
+
+  else if(scenario ==2)
+    {
+      phiDir = phithetamin;
+    }
+  else 
+    return kFALSE;
+
+
+      Phi = phiTA;
+
+      if( TMath::Abs(CalcdPhi(phiDir, phiTA)) <TMath::Pi()/2)
+       return kTRUE;
+      else {
+       Phi+=TMath::Pi();
+       if(Phi>TwoPi) Phi-=TwoPi;
+       return kTRUE;
+      }
+
+
+ return kTRUE;
+} 
+*/
+
+Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxes(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
+{
+  // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
+
+  Double_t xmean , ymean;
+
+  Double_t TwoPi = TMath::TwoPi();
+
+  xmean = 0;
+  ymean = 0;
+  Phi= -999;
+  //find axes
+
+  Double_t x=0;
+  Double_t x2=0;
+  Double_t y=0; 
+  Double_t y2=0;
+  Double_t xy=0;
+  
+  //  Double_t phiLoc, thetaLoc;
+
+  Int_t N=0;
+
+  Int_t Nev =  list.GetSize();
+  if( Nev <2) return kFALSE;
+
+
+
+
+  Int_t Index = -1;
+  Double_t val = -1;
+  if(scenario==2) val = 100;
+  Double_t phiDir = -99;
+
+
+
+  for(Int_t ip=0; ip< Nev; ip++) {
+    if(list.At(ip)<0) break;
+
+    TVector3* part = (TVector3*) GetAt(list.At(ip));
+    if(!part) continue;
+
+
+    Double_t xx, yy;
+    /*
+    TVector3 vtmp(*part);
+    GetLocalMom(vec, &vtmp);
+    Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
+    Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
+    */
+
+    xx = CalcdPhi(part->Phi(), vec.Phi());
+    yy = part->Eta() - vec.Eta();
+
+    x+=xx;
+    y+=yy;
+    x2+=(xx*xx);
+    y2+=(yy*yy);
+    xy+=(xx*yy);
+    N++;
+
+    /*
+  if(scenario ==0)
+    {
+      if(vtmp.Mag()>val)
+       {
+         val=vtmp.Mag();
+         Index = ip;
+         phiDir=vtmp.Phi();
+       }
+    }
+  else if(scenario ==1)
+    {
+      if(vtmp.Z()>val)
+       {
+         val=vtmp.Z();
+         Index = ip;
+         phiDir=vtmp.Phi();
+       }
+    }
+  else if(scenario ==2)
+    {
+      if(vtmp.Theta()<val)
+       {
+         val=vtmp.Theta();
+         Index = ip;
+         phiDir=vtmp.Phi();
+       }
+    }
+  else 
+*/
+if(scenario ==3)
+    {
+      if(part->Pt()> val)
+       {
+         val =part->Pt();
+         Index = ip;
+         //      phiDir=vtmp.Phi();
+         phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
+       }
+    }
+
+
+  }
+
+  if(N<2) return kFALSE;
+
+  //  return kFALSE;
+
+  x/=N;
+  y/=N;
+  x2/=N;
+  y2/=N;
+  xy/=N;
+
+  
+  Double_t A = xy-x*y;
+  Double_t B = x2-x*x+y*y-y2;
+  Double_t D = TMath::Sqrt(B*B+4*A*A);
+
+  //  printf("N = %d  A = %f\n",N, A);
+  //  list.Print();
+
+  Double_t a1 = (-B+D)/2/A;
+  //  Double_t a2 = (-B-D)/2/A;
+  //  Double_t b1 = y-a1*x;
+  //  Double_t b2 = y-a2*x;
+  //  Double_t phiDir = TMath::ATan2(y, x);
+  // while(phiDir<0) phiDir+=TwoPi;
+    Double_t phi = TMath::ATan(a1);
+    
+  while(phi>TwoPi) phi-=TwoPi;
+  while(phi< 0 )   phi+=TwoPi;
+
+  Phi=CalcdPhi0To2pi(phi, 0);
+
+  if(scenario!=4 && Index<0) return kFALSE;
+
+  if(Index>=0)
+    {
+      if( TMath::Abs(CalcdPhi(phiDir, phi)) <TMath::Pi()/2)
+       return kTRUE;
+      else
+       {
+         phi+=TMath::Pi();
+         Phi = CalcdPhi0To2pi(phi);
+         return kTRUE;
+       }
+    }
+
+
+  Double_t xmax = -100;
+  Double_t xmin =  100;
+
+  xmean = x;
+  ymean = y;
+
+
+  for(Int_t ip=0; ip< Nev; ip++) {
+    if(list.At(ip)<0) break;
+
+    TVector3* part = (TVector3*) GetAt(list.At(ip));
+    if(!part) continue;
+
+    TVector3 vtmp(*part);
+    GetLocalMom(vec, &vtmp);
+
+    Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
+    Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
+
+    Double_t x1 = (xx-xmean)*TMath::Cos(Phi) +(yy-ymean)*TMath::Sin(Phi);
+    if(x1 > xmax) xmax = x1;
+    if(x1 < xmin) xmin = x1;
+    //    printf("c3\n");
+  }
+   
+  if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
+
+
+
+   while(Phi>TwoPi) Phi-=TwoPi;
+
+  return kTRUE;
+}
+
+
+
+
+Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxesCorr(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario, Double_t R)
+{
+  // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
+
+  Double_t xmean , ymean;
+
+  Double_t TwoPi = TMath::TwoPi();
+
+  xmean = 0;
+  ymean = 0;
+  Phi= -999;
+  //find axes
+
+  Double_t x=0;
+  Double_t x2=0;
+  Double_t y=0; 
+  Double_t y2=0;
+  Double_t xy=0;
+  
+  //  Double_t phiLoc, thetaLoc;
+
+  Int_t N=0;
+
+  Int_t Nev =  list.GetSize();
+  if( Nev <2) return kFALSE;
+
+
+
+
+  Int_t Index = -1;
+  Double_t val = -1;
+  if(scenario==2) val = 100;
+  Double_t phiDir = -99;
+
+
+
+  for(Int_t ip=0; ip< Nev; ip++) {
+    if(list.At(ip)<0) break;
+
+    TVector3* part = (TVector3*) GetAt(list.At(ip));
+    if(!part) continue;
+
+    Double_t xx, yy;
+    /*
+    TVector3 vtmp(*part);
+    GetLocalMom(vec, &vtmp);
+    Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
+    Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
+    */
+
+    xx = CalcdPhi(part->Phi(), vec.Phi());
+    yy = part->Eta() - vec.Eta();
+
+    x+=xx;
+    y+=yy;
+    x2+=(xx*xx);
+    y2+=(yy*yy);
+    xy+=(xx*yy);
+    N++;
+
+    /*
+  if(scenario ==0)
+    {
+      if(vtmp.Mag()>val)
+       {
+         val=vtmp.Mag();
+         Index = ip;
+         phiDir=vtmp.Phi();
+       }
+    }
+  else if(scenario ==1)
+    {
+      if(vtmp.Z()>val)
+       {
+         val=vtmp.Z();
+         Index = ip;
+         phiDir=vtmp.Phi();
+       }
+    }
+  else if(scenario ==2)
+    {
+      if(vtmp.Theta()<val)
+       {
+         val=vtmp.Theta();
+         Index = ip;
+         phiDir=vtmp.Phi();
+       }
+    }
+  else 
+*/
+if(scenario ==3)
+    {
+      if(part->Pt()> val)
+       {
+         val =part->Pt();
+         Index = ip;
+         //      phiDir=vtmp.Phi();
+         phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
+       }
+    }
+
+
+  }
+
+  if(N<2) return kFALSE;
+  if(scenario!=4 && Index<0) return kFALSE;
+
+  //  return kFALSE;
+
+  x/=N;
+  y/=N;
+  x2/=N;
+  y2/=N;
+  xy/=N;
+
+  
+  Double_t A = xy-x*y;
+  Double_t B = x2-x*x+y*y-y2;
+  Double_t D = TMath::Sqrt(B*B+4*A*A);
+
+  //  printf("N = %d  A = %f\n",N, A);
+  //  list.Print();
+
+  Double_t a1 = (-B+D)/2/A;
+  //  Double_t a2 = (-B-D)/2/A;
+  //  Double_t b1 = y-a1*x;
+  //  Double_t b2 = y-a2*x;
+  //  Double_t phiDir = TMath::ATan2(y, x);
+  // while(phiDir<0) phiDir+=TwoPi;
+    Double_t phi = TMath::ATan(a1);
+    
+  Phi=CalcdPhi0To2pi(phi, 0);
+  if( TMath::Abs(CalcdPhi(phiDir, phi))  > TMath::Pi()/2)
+    Phi=CalcdPhi0To2pi(phi+TMath::Pi(), 0);
+
+
+
+
+  if(Index>=0)
+    {
+      if(N==2) 
+       return kTRUE;
+
+
+
+      Double_t phiMinus = 0;
+      Int_t Nminus = 0;
+      for(Int_t ip=0; ip< Nev; ip++) 
+       {
+         if(list.At(ip)<0) break;
+
+         TVector3* part = (TVector3*) GetAt(list.At(ip));
+         if(!part) continue;
+         Double_t xx0 = CalcdPhi(part->Phi(), vec.Phi());
+         Double_t yy0 = part->Eta() - vec.Eta();
+
+         Double_t xx  = (x*N - xx0)/(N-1);
+         Double_t yy  = (y*N - yy0)/(N-1);
+         Double_t xx2 = (x2*N - xx0*xx0)/(N-1);
+         Double_t yy2 = (y2*N - yy0*yy0)/(N-1);
+         Double_t xxyy= (xy*N - xx0*yy0)/(N-1);
+
+
+         Double_t AA = xxyy-xx*yy;
+         Double_t BB = xx2-xx*xx+yy*yy-yy2;
+         Double_t DD = TMath::Sqrt(BB*BB+4*AA*AA);
+         Double_t aa1 = (-BB+DD)/2/AA;
+         Double_t phi1 = TMath::ATan(aa1);
+    
+         if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
+           phi1+=TMath::Pi();
+
+         phiMinus+=CalcdPhi(Phi, phi1);
+         Nminus++;
+
+       }
+
+
+      Double_t phiPlus = 0;
+      Int_t Nplus = 0;
+      for(Int_t ip=0; ip< -Nev; ip++) 
+       {
+         if(list.At(ip)<0) break;
+
+         TVector3* part = (TVector3*) GetAt(list.At(ip));
+         if(!part) continue;
+         Double_t rtmp = gRandom->Uniform(0, R);
+         Double_t phitmp = gRandom->Uniform(0, TMath::TwoPi());
+         Double_t xx0 = rtmp*TMath::Cos(phitmp);
+         Double_t yy0 = rtmp*TMath::Sin(phitmp);
+
+         Double_t xx  = (x*N + xx0)/(N+1);
+         Double_t yy  = (y*N + yy0)/(N+1);
+         Double_t xx2 = (x2*N + xx0*xx0)/(N+1);
+         Double_t yy2 = (y2*N + yy0*yy0)/(N+1);
+         Double_t xxyy= (xy*N + xx0*yy0)/(N+1);
+
+
+         Double_t AA = xxyy-xx*yy;
+         Double_t BB = xx2-xx*xx+yy*yy-yy2;
+         Double_t DD = TMath::Sqrt(BB*BB+4*AA*AA);
+         Double_t aa1 = (-BB+DD)/2/AA;
+         Double_t phi1 = TMath::ATan(aa1);
+    
+         if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
+           phi1+=TMath::Pi();
+
+         phiPlus+=CalcdPhi(Phi, phi1);
+         Nplus++;
+
+       }
+
+
+      Phi+=((phiMinus+phiPlus)/(Nminus+Nplus+1));
+
+      if( TMath::Abs(CalcdPhi(phiDir, Phi)) > TMath::Pi()/2)
+       Phi+=TMath::Pi();
+
+      Phi=CalcdPhi0To2pi(Phi, 0);
+
+      return kTRUE;
+    }
+
+
+  Double_t xmax = -100;
+  Double_t xmin =  100;
+
+  xmean = x;
+  ymean = y;
+
+
+  for(Int_t ip=0; ip< Nev; ip++) {
+    if(list.At(ip)<0) break;
+
+    TVector3* part = (TVector3*) GetAt(list.At(ip));
+    if(!part) continue;
+
+    TVector3 vtmp(*part);
+    GetLocalMom(vec, &vtmp);
+
+    Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
+    Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
+
+    Double_t x1 = (xx-xmean)*TMath::Cos(Phi) +(yy-ymean)*TMath::Sin(Phi);
+    if(x1 > xmax) xmax = x1;
+    if(x1 < xmin) xmin = x1;
+    //    printf("c3\n");
+  }
+   
+  if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
+
+
+
+   while(Phi>TwoPi) Phi-=TwoPi;
+
+  return kTRUE;
+}
+
+
+/////////////////////////////////////
+
+ TH2F *fhPhiEtaTrack;//! 
+
+ TH1F *fhTMA_JAA[3];//! 
+ TH1F *fhTMA_JAp[3];//! 
+ TH1F *fhTMA_B1AA[3];//! 
+ TH1F *fhTMA_B1Ap[3];//! 
+ TH1F *fhTMA_B2AA[3];//! 
+ TH1F *fhTMA_B2Ap[3];//! 
+
+
+ Int_t     fPtJetNbin;
+ TArrayD   fPtJetArray;
+ Int_t     fPsiVsRNbin;
+ Double_t  fRmax;
+ Int_t     fPsiVsPhiNbin;
+
+ UInt_t   fFilterMask;
+ Double_t fEtaTrackMax;
+ Double_t fPtTrackMin;
+ Double_t fPtTrackMax;
+ Double_t fPtJmin;
+ Double_t fPtJmax;
+ Double_t fPtJ;
+
+ TVector3 fJvec;
+
+
+  AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AliAnalysisTaskJetShapeHM(TString comment, Bool_t kMC):TObject(),
+fComment(comment),
+fkMC(kMC),
+fhEvent(0),
+fhMult(0),
+fhPtJ(0),
+fhPtJvsPtCorr(0),
+fhPsiVsR(0),
+fhPsiVsRPtJ(0),
+fhPhiEtaTrack(0),
+fPtJetNbin(0),
+fPtJetArray(0),
+fPsiVsRNbin(0),
+fRmax(0),
+fPsiVsPhiNbin(0),
+fFilterMask(0),
+fEtaTrackMax(0),
+fPtTrackMin(0),
+fPtTrackMax(0),
+fPtJmin(0),
+fPtJmax(0),
+fPtJ(0),
+fJvec(0,0,0)
+{
+
+  for(Int_t i=0; i<3; i++)
+    {
+      fhTMA_JAA[i]=0;
+      fhTMA_JAp[i]=0;
+      fhTMA_B1AA[i]=0;
+      fhTMA_B1Ap[i]=0;
+      fhTMA_B2AA[i]=0;
+      fhTMA_B2Ap[i]=0;
+    }
+
+
+  fPtJetNbin = 0;
+  fPtJetArray = 0;
+
+  SetRBins();
+  SetPhiNbins();
+  SetFilterMask();
+  SetEtaTrackMax();
+  SetPtTrackRange();
+  SetPtJetRange();
+
+  fJvec.SetXYZ(0,0,0);
+  fPtJ = 0;
+}
+
+
+AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::~AliAnalysisTaskJetShapeHM()
+{
+  /*
+  if(fPtJetArray)
+    delete [] fPtJetArray;
+  fPtJetArray=0;
+  */
+}
+
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::SetPtJetBins(Int_t Nbin, Double_t *array)
+{
+
+  /*
+  if(fPtJetArray)
+    delete [] fPtJetArray;
+  fPtJetArray=0;
+
+  fPtJetNbin = Nbin;
+  fPtJetArray = new Double_t[fPtJetNbin+1];
+  for(Int_t i=0; i<= fPtJetNbin; i++) fPtJetArray[i]= array[i]; 
+  */
+
+
+  fPtJetNbin = Nbin;
+  fPtJetArray.Set(fPtJetNbin+1, array);
+}
+
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::InitHistos()
+{
+
+  fhEvent    = Hist1D("hEvent" , 3  , -0.5,   2.5, "Event");
+  fhMult     = Hist1D("hMult"  , 101, -0.5, 100.5, "N_{ch}");
+
+  fhPhiEtaTrack = Hist2D("hPhiEtaTrack", 100, 0, TMath::TwoPi(), 20, -1, 1., "#phi_{tr}", "#eta_{tr}");
+
+  if(fPtJetNbin<1) {
+    Int_t Nbin = 13;
+    Double_t array[]= {0., 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 130, 160, 200};
+    SetPtJetBins(Nbin, array);
+  }
+
+  fhPtJ      = Hist1D("hPtJ"  , fPtJetNbin, fPtJetArray.GetArray(), "Pt_{J} GeV/c");
+  fhPsiVsR   = Hist1D("hPsiVsR", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
+  fhPsiVsRPtJ   = Hist2D("hPsiVsRPtJ", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
+
+  fhPtJvsPtCorr  = Hist2D("fhPtJvsPtCorr", fPtJetNbin, fPtJetArray.GetArray(), 100, -100, 100, "P_{tJ} GeV/c", "P_{tJ} - P_{tJ,corr} GeV/c" , 1);
+
+  for(Int_t i=0; i<3; i++)
+    {
+      fhTMA_JAA[i]  = Hist1D(Form("fhTMA_JAA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
+      fhTMA_JAp[i]  = Hist1D(Form("fhTMA_JAp_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
+
+      fhTMA_B1AA[i]  = Hist1D(Form("fhTMA_B1AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
+      fhTMA_B1Ap[i]  = Hist1D(Form("fhTMA_B1Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
+
+      fhTMA_B2AA[i]  = Hist1D(Form("fhTMA_B2AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
+      fhTMA_B2Ap[i]  = Hist1D(Form("fhTMA_B2Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
+    }
+
+
+}
+
+
+
+
+
+void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
+{
+  list->Add(fhEvent);
+  list->Add(fhMult);
+  list->Add(fhPhiEtaTrack);
+  list->Add(fhPtJ);
+  list->Add(fhPsiVsR);
+  list->Add(fhPsiVsRPtJ);
+  list->Add(fhPtJvsPtCorr);
+
+  for(Int_t i=0; i<3; i++)
+    {
+      list->Add(fhTMA_JAA[i]);
+      list->Add(fhTMA_JAp[i]);
+
+      list->Add(fhTMA_B1AA[i]);
+      list->Add(fhTMA_B1Ap[i]);
+
+      list->Add(fhTMA_B2AA[i]);
+      list->Add(fhTMA_B2Ap[i]);
+    }
+
+
+}
+
+
+
+
+Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE,  AliAODJet *jet, Double_t DeltaPtJ)
+{
+
+
+
+  Int_t size = aodE->GetNumberOfTracks();
+
+  TClonesArray *arrP = 0;
+
+  if(fkMC)
+    {
+      arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
+     if(!arrP) 
+       {
+        printf("ERROR: no Info about particles in AOD\n");
+        return kFALSE;
+       }
+     size = arrP->GetEntriesFast();
+    }
+
+
+  Int_t IndexArray[size];
+
+  TClonesArray farray("TVector3", size);
+
+  Int_t counter = 0;
+  Int_t counterAll = 0;
+  Double_t pJ[3] = {0, 0, 0};
+
+
+  TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
+
+  for(int it = 0;it < size;++it){
+
+    TVector3 vec;
+
+    if(fkMC)
+      {
+       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
+       if(!part)continue;
+       if(!part->IsPhysicalPrimary())continue;
+       if(part->Charge()==0)continue;
+       vec.SetXYZ(part->Px(), part->Py(), part->Pz());
+      }
+    else 
+      {
+       AliAODTrack *tr = aodE->GetTrack(it);
+       if(!tr) continue;
+       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+       vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
+      }
+
+
+    if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
+    if(vec.Pt()< fPtTrackMin)continue;
+    if(vec.Pt()> fPtTrackMax) {return kFALSE;}
+
+
+    new(farray[counterAll]) TVector3(vec);
+    counterAll++;
+
+
+    Double_t R = CalcR(vecJ, vec);
+    if(R> fRmax) continue;
+
+    pJ[0]+=vec.Px();
+    pJ[1]+=vec.Py();
+    pJ[2]+=vec.Pz();
+
+    IndexArray[counter] = it;
+    counter++;
+  }
+
+  fhMult->Fill(counter);
+  if(counter<1) return kFALSE;
+
+  fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
+
+  fPtJ =  TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
+  if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
+  fhPtJ->Fill(fPtJ);
+
+  fhPtJvsPtCorr->Fill(fPtJ, jet->Px() - DeltaPtJ);
+
+  for(Int_t it = 0; it<counter; it++)
+    {
+      TVector3 vec;
+
+      if(fkMC)
+       {
+         AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
+         vec.SetXYZ(part->Px(), part->Py(), part->Pz());
+       }
+      else 
+       {
+         AliAODTrack *tr = aodE->GetTrack(IndexArray[it]);
+         vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
+       }
+
+       Double_t R = CalcR(fJvec, vec);
+       //      Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
+       Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
+       fhPsiVsR->Fill(R, probL);
+       fhPsiVsRPtJ->Fill(R, fPtJ, probL);
+       Double_t phi = vec.Phi();
+       while(phi<0) phi+=TMath::TwoPi();
+       while(phi>TMath::TwoPi()) phi-=TMath::TwoPi();
+       fhPhiEtaTrack->Fill(phi, vec.Eta());
+
+     }
+
+
+  // ang. distr.
+  //     fMyTool->Clean();
+  AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
+
+      fMyTool->SetNEntries(counterAll);
+      fMyTool->SetVecJ(vecJ);
+      fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
+      fMyTool->Init();
+      Int_t scenario = 3;
+
+      /*
+       to be added!!!!!!!!!
+      fhTMA_B1AA[i]=0;
+      fhTMA_B2AA[i]=0;
+      */
+
+      for(Int_t l=0; l<3; l++)
+       {
+         Double_t phiA, phi1;
+
+         //      Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
+         //      Double_t ptRJ  = fMyTool->GetPRecJ(l,0).Pt();
+         Int_t N1 = fMyTool->GetSizeJ(l,1);
+         Double_t dphi = -999;
+
+
+         if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0)  , vecJ, phiA, scenario))
+           continue;
+
+         //      f2JEvent->SetPhiJL(l,0, phiA);
+
+         if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1)  , vecJ, phi1, scenario))
+           {
+             fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
+           }
+
+
+         for(Int_t ip =0; ip<N1; ip++) 
+           {
+             phi1 = fMyTool->GetLocPhiJ(l,1,ip);
+             dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
+             fhTMA_JAp[l]->Fill(dphi);
+           }
+
+
+         Int_t NB1 = fMyTool->GetSizeB1(l, 1);
+         for(Int_t ip =0; ip<NB1; ip++) 
+           {
+             phi1 = fMyTool->GetLocPhiB1(l,1,ip);
+             dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
+             fhTMA_B1Ap[l]->Fill(dphi);
+           }
+
+         Int_t NB2 = fMyTool->GetSizeB2(l, 1);
+         for(Int_t ip =0; ip<NB2; ip++) 
+           {
+             phi1 = fMyTool->GetLocPhiB2(l,1,ip);
+             dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
+             fhTMA_B2Ap[l]->Fill(dphi);
+           }
+
+       }
+
+    fhEvent->Fill(1);
+
+    delete fMyTool;
+
+  return kTRUE;
+}
+
+
+Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
+{
+
+  Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
+  Double_t deta = v1.Eta() - v2.Eta();
+  Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
+  return RB;
+}
+
+
+Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
+{
+
+  while(phi1<0) phi1+=TMath::TwoPi();
+  while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
+
+  while(phi2<0) phi2+=TMath::TwoPi();
+  while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
+
+  Double_t dphi = phi1- phi2;
+  if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
+  if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
+
+  return dphi;
+}
+
+
+
+
+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)
+{
+// create a 1D histogram
+
+  TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (yLabel) res->GetYaxis()->SetTitle(yLabel);
+  res->SetLineColor(color);
+  return res;
+}
+
+
+TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray,  const char* xLabel, Int_t color, const char* yLabel)
+{
+// create a 1D histogram
+
+  TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (yLabel) res->GetYaxis()->SetTitle(yLabel);
+  res->SetLineColor(color);
+  return res;
+}
+
+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)
+{
+// create a 2D histogram
+
+  TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (yLabel) res->GetYaxis()->SetTitle(yLabel);
+  //  res->SetMarkerStyle(kFullCircle);
+  //  res->SetOption("E");
+  res->SetLineColor(color);
+  //  fOutputList->Add(res);
+  return res;
+}
+
+
+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)
+{
+// create a 2D histogram
+
+  TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (yLabel) res->GetYaxis()->SetTitle(yLabel);
+  if (zLabel) res->GetZaxis()->SetTitle(zLabel);
+  //  res->SetMarkerStyle(kFullCircle);
+  //  res->SetOption("E");
+  res->SetLineColor(color);
+  //  fOutputList->Add(res);
+  return res;
+}
+
+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)
+{
+// create a 2D histogram
+
+  TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (yLabel) res->GetYaxis()->SetTitle(yLabel);
+  if (zLabel) res->GetZaxis()->SetTitle(zLabel);
+  //  res->SetMarkerStyle(kFullCircle);
+  //  res->SetOption("E");
+  res->SetLineColor(color);
+  //  fOutputList->Add(res);
+  return res;
+}
+
+
+
+//////////////////////////////////////
+
+
+/*
+//________________________________________________________________________
+void AnalysisJetMain::ConnectInputData(Option_t *)
+{
+  // Connect ESD
+  // Called once
+
+  fESD=dynamic_cast<AliESDEvent*>(InputEvent());
+  //   if (!fESD) {
+  //     AliError("ESD not available");
+  //     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
+  //     fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
+  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
+  //  assume that the AOD is in the general output...
+  fAODOut  = AODEvent();
+
+
+}
+*/
+
+void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
+{
+   fJetBranchName[0] = branch1;
+   fJetBranchName[1] = branch2;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskJetShape::UserCreateOutputObjects()
+{
+  //printf("Open1\n");
+  //       const char *nameF = OpenFile(1)->GetName();
+  //printf("Open2 %s\n",nameF);
+
+  //  fTriggerAnalysis = new AliTriggerAnalysis();
+  //  if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
+
+
+  fOutputList = new TList();
+  fOutputList->SetOwner();
+
+  fhPtJL  = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
+  fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
+
+  
+  fanalJetSubStrHM->SetFilterMask(fFilterMask);
+  fanalJetSubStrHM->InitHistos();
+  fanalJetSubStrHM->AddToList(fOutputList);
+
+  printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
+
+  if(fkMC)
+    {
+
+      printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
+      fanalJetSubStrMCHM->InitHistos();
+      fanalJetSubStrMCHM->AddToList(fOutputList);
+    }
+  
+
+
+
+  PostData(1, fOutputList);
+
+  /*
+   OpenFile(1);
+
+  fOutputTree = new TTree("tree","Tree with Ali2JEvent");
+  fOutputTree->Branch("event",&f2JEvent);
+  */
+
+//    PostData(1, fOutputTree);
+}
+
+/*
+Bool_t AliAnalysisTaskJetSpectrum2::Notify()
+{
+
+  // Fetch the aod also from the input in,
+  // have todo it in notify
+  
+  
+  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
+  //  assume that the AOD is in the general output...
+  fAODOut  = AODEvent();
+  
+  if(fNonStdFile.Length()!=0){
+    // case that we have an AOD extension we need can fetch the jets from the extended output
+    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+    fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
+    if(!fAODExtension){
+      if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
+    }
+  }
+
+
+
+  return kTRUE;
+}
+*/
+
+
+
+
+//________________________________________________________________________
+void AliAnalysisTaskJetShape::UserExec(Option_t *) 
+{
+  //   return;
+//    if(f2JEvent)
+//      delete f2JEvent;
+//    f2JEvent = new Ali2JEvent();
+//    return;
+  
+  //  return;
+  if(fDebug) printf("NEW EVENT 0\n");
+
+  if(!IsGoodEvent()) return;
+
+  if(fDebug) printf("NEW EVENT 1\n");
+
+  //   printf("\n\n\n NEW EVENT\n");
+ if(fDebug)  Printf("Event #%5d", (Int_t) fEntry);
+
+   AliAODEvent* aodE = 0;
+   if(!aodE){
+     if(!fESD)aodE = fAODIn;
+     else aodE = fAODOut;}
+
+
+
+  // centrality selection
+   AliCentrality *cent = 0x0;
+   Double_t centrality = 0.; 
+
+   if(fESD) {cent = fESD->GetCentrality();
+     if(cent) centrality = cent->GetCentralityPercentile("V0M");}
+   else     centrality=aodE->GetHeader()->GetCentrality();
+
+
+   if(!fkIsPbPb) {
+     centrality=1.;
+   }
+
+  if(fDebug) printf("NEW EVENT 2\n");
+
+   Bool_t  fUseAOD049 = kFALSE;// kTRUE;// 
+      if(fUseAOD049&&centrality>=0){
+       Float_t v0=0;
+       AliAODVZERO* aodV0 = aodE->GetVZEROData();
+       v0+=aodV0->GetMTotV0A();
+       v0+=aodV0->GetMTotV0C();
+       if(centrality==0 && v0 < 19500) return ;//filtering issue
+       Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
+       Float_t val= 1.30552 +  0.147931 * v0;
+       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, 
+            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, 
+            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, 
+            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, 
+            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, 
+            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, 
+            13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
+       if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] ) 
+         return; //outlier
+      }
+   
+
+    if(fDebug) printf("centrality: %f\n", centrality);
+    if (centrality < fCentMin || centrality > fCentMax){
+    //    PostData(1, fOutputList);
+    return;
+    }
+
+   //   fhNEvent->Fill(0); 
+
+   // accepted events  
+   // -- end event selection --
+  
+    //  aodE->Print();
+    //   aodE->GetList()->Print();
+   
+   // get background
+   AliAODJetEventBackground* externalBackground = 0;
+   AliAODJetEventBackground* externalBackgroundMC = 0;
+
+
+   Float_t rho = 0;
+   Float_t rho_MC=0;
+
+   if(fkIsPbPb)
+     {
+       if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
+        externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
+        if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
+       }
+       if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
+        externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
+        if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
+       }
+       if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
+        externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
+        if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
+       }
+       if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
+        externalBackgroundMC =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
+        if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
+       }
+       if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
+        externalBackgroundMC =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
+        if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
+       }
+       if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
+        externalBackgroundMC =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
+        if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
+       }
+       if(externalBackground)rho = externalBackground->GetBackground(0);
+       if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
+     }
+
+
+
+
+   //   printf("rho = %f\n",rho);
+
+    //   return;
+
+   TClonesArray *Jets[2];
+   Jets[0]=0;
+   Jets[1]=0;
+   if(fAODOut&&!Jets[0]){
+   Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); 
+   Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
+   if(fAODExtension && !Jets[0]){ 
+   Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
+   Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
+   if(fAODIn&&!Jets[0]){
+   Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data())); 
+   Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
+
+
+   if(!Jets[0]) {
+     Printf("no friend!!!\n");
+     return;
+   }
+  Int_t nJ = Jets[0]->GetEntries();
+  Float_t ptmax = 0.;
+  Int_t   imax  = -1;
+
+  if(fDebug) printf("NEW EVENT 3\n");
+
+
+//
+// Find highest pT jet with pt > 20 GeV
+//
+//  fPtJcorrMin=30;
+
+  Float_t etaJmax = 0.4;
+
+     for (Int_t i = 0; i < nJ; i++) {
+       AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
+       if (!jet) continue;
+       //       jet->Print("0");
+       Float_t ptJ  = jet->Pt();
+       Float_t etaJ = TMath::Abs(jet->Eta());
+
+       Float_t area = jet->EffectiveAreaCharged();
+       Float_t ptJcorr=ptJ-rho*area;
+       
+       if ((ptJcorr > fPtJcorrMin) && (ptJcorr  > ptmax) && etaJ < etaJmax) {
+        ptmax = ptJcorr;
+        imax = i;
+       }
+       
+     }
+
+     TVector3 vecJ;
+
+     AliAODJet* jetL = 0;
+
+     Float_t areaJL, ptJLcorr;
+
+     if (imax != -1)  {
+
+     jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
+     areaJL   = jetL->EffectiveAreaCharged();
+     ptJLcorr  = jetL->Pt()-rho*areaJL;
+
+     fhPtJL->Fill(ptJLcorr);
+     fhAreaJL->Fill(areaJL);
+     vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
+     fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
+     }
+
+
+
+     if(!fkMC)
+       {
+        PostData(1, fOutputList);
+        return;
+       }
+
+
+     nJ = Jets[1]->GetEntries();
+     ptmax = 0;
+     imax = -1;
+     Float_t areaJL_MC=0;
+
+     for (Int_t i = 0; i < nJ; i++) {
+       AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
+       if (!jet) continue;
+       Float_t ptJ1  = jet->Pt();
+       Float_t etaJ1 = TMath::Abs(jet->Eta());
+
+       areaJL_MC = jet->EffectiveAreaCharged();
+       Float_t ptJcorr=ptJ1-rho*areaJL_MC;
+
+       //    jet->Print();
+
+       if ((ptJcorr > fPtJcorrMin) && (ptJcorr  > ptmax) && etaJ1 < etaJmax) {
+        ptmax = ptJcorr;
+        imax = i;
+       }
+       
+     }
+
+     AliAODJet* jetMCL=0;
+
+    if (imax != -1)  {
+      jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
+      fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
+     }
+
+
+
+
+
+     PostData(1, fOutputList);
+     return;
+
+
+
+
+}      
+
+Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
+{
+  while(phi1<0) phi1+=TMath::TwoPi();
+  while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
+
+  while(phi2<0) phi2+=TMath::TwoPi();
+  while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
+
+  Double_t dphi = phi1- phi2;
+  if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
+  if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
+
+  //  Double_t dphi = phi1- phi2;
+  //  while(dphi<0) dphi+=TMath::TwoPi();
+  //  while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
+  return dphi;
+}
+
+Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
+{
+
+  Double_t dphi = CalcdPhi(phi1, phi2);
+  while(dphi<0) dphi+=TMath::TwoPi();
+  while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
+  return dphi;
+}
+
+
+Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
+{
+  
+  fESD=dynamic_cast<AliESDEvent*>(InputEvent());
+  //   if (!fESD) {
+  //     AliError("ESD not available");
+  //     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
+  //     fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
+  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
+  //  assume that the AOD is in the general output...
+  fAODOut  = AODEvent();
+  
+    
+   static AliAODEvent* aod = 0;
+   // take all other information from the aod we take the tracks from
+   if(!aod){
+     if(!fESD)aod = fAODIn;
+     else aod = fAODOut;}
+     
+
+   if(fNonStdFile.Length()!=0){
+    // case that we have an AOD extension we need can fetch the jets from the extended output
+     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
+   }
+
+   if(fDebug){
+    if(fAODIn) printf("fAODIn\n");
+    if(fAODOut) printf("fAODOut\n");
+    if(fAODExtension) printf("fAODExtension\n");
+   }
+
+   // -- event selection --
+
+   // physics selection
+   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+   ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+   //   cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
+   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
+      if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
+      return kFALSE;
+   }
+
+   // vertex selection
+   //   if(!aod){
+   //     if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
+   //     fhNEvent->Fill(3);
+   //      PostData(1, fOutputList);
+   //     return kFALSE;
+   //   }
+
+   AliAODVertex* primVtx = aod->GetPrimaryVertex();
+
+   if(!primVtx){
+     if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
+     return kFALSE;
+   }
+
+   Int_t nTracksPrim = primVtx->GetNContributors();
+   if ((nTracksPrim < fVtxMinContrib) ||
+         (primVtx->GetZ() < fVtxZMin) ||
+         (primVtx->GetZ() > fVtxZMax) ){
+      if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
+      return  kFALSE;
+   }
+
+   /*
+   // event class selection (from jet helper task)
+   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
+   if(fDebug) Printf("Event class %d", eventClass);
+   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
+      return  kFALSE;
+   }
+   */
+    return kTRUE;
+}
+
+
+
+
+
+
+//________________________________________________________________________
+void AliAnalysisTaskJetShape::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+ //    fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
+
+  printf("MyTaskTestTrigger::Terminate()\n\n\n");
+
+  
+  //  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  //  fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
+
+}
+
+
+
+/*
+Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd,  Double_t Vxyz[3], Int_t type)
+{
+  // type
+  //     0 -> SPD vtx
+  //     1 -> ESD vtx
+  const AliESDVertex* vtx = 0;
+
+  if(type==0) {
+    vtx = esd->GetPrimaryVertexSPD();
+    if(!vtx) return kFALSE;
+    if(vtx->GetNContributors() < 1) return kFALSE;
+    TString vtxTyp = vtx->GetTitle();
+    if (vtxTyp.Contains("vertexer: Z")) {
+      if (vtx->GetDispersion()>0.04) return kFALSE;
+      if (vtx->GetZRes()>0.25) return kFALSE;
+    }
+  }
+  if(type==1) {
+    vtx = esd->GetPrimaryVertexTracks();
+    if(!vtx) return kFALSE;
+    if(vtx->GetNContributors() < 1) return kFALSE;
+  }
+
+   Vxyz[0] = vtx->GetXv();
+   Vxyz[1] = vtx->GetYv();
+   Vxyz[2] = vtx->GetZv();
+   return kTRUE;
+}
+*/
+
+
+
+
+
+
+
+
+
+TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax,  const char* xLabel, Int_t color)
+{
+// create a 1D histogram
+
+  TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
+  if (xLabel) h->GetXaxis()->SetTitle(xLabel);
+  //  res->SetMarkerStyle(kFullCircle);
+  //  res->SetOption("E");
+  h->SetLineColor(color);
+  //  fOutputList->Add(h);
+  return h;
+}
+
+
+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)
+{
+// create a 2D histogram
+
+  TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (xLabel) res->GetYaxis()->SetTitle(yLabel);
+  //  res->SetMarkerStyle(kFullCircle);
+  //  res->SetOption("E");
+  res->SetLineColor(color);
+  //  fOutputList->Add(res);
+  return res;
+}
+
+
+TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, 
+                                                  Int_t nBinsy, Double_t yMin, Double_t yMax, 
+                                                 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
+{
+// create a 3D histogram
+
+  TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
+  if (xLabel) res->GetXaxis()->SetTitle(xLabel);
+  if (yLabel) res->GetYaxis()->SetTitle(yLabel);
+  if (zLabel) res->GetZaxis()->SetTitle(zLabel);
+  //  res->SetMarkerStyle(kFullCircle);
+  //  res->SetOption("E");
+  res->SetLineColor(color);
+  //  fOutputList->Add(res);
+  return res;
+}
+
+
+
+
+//#endif
+
+
+
+
+
+
+
+
diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetShape.h b/PWGJE/UserTasks/AliAnalysisTaskJetShape.h
new file mode 100644 (file)
index 0000000..8055b61
--- /dev/null
@@ -0,0 +1,404 @@
+#ifndef ALIANALYSISTASKJETSHAPE_h
+#define ALIANALYSISTASKJETSHAPE_h
+
+#include "AliAnalysisTaskSE.h"
+#include "AliVEvent.h"
+//#include "AliTriggerAnalysis.h"
+//#include "THnSparse.h"
+#include "TMath.h"
+#include "TString.h"
+#include "TH1F.h"
+#include "TClonesArray.h"
+
+class TH2F;
+class TH3F;
+class TVector3;
+
+class AliAODExtension;
+class TNtuple;
+class TTree;
+class TList;
+class AliAODEvent;
+class AliAODJet;
+class TArrayD;
+
+
+
+
+
+
+
+class AliAnalysisTaskJetShape : public AliAnalysisTaskSE {
+
+ public:
+
+
+class AliAnalysisTaskJetShapeTool :public TObject {
+public:
+AliAnalysisTaskJetShapeTool();
+//AnalJetSubStrTool(TList *list, TVector3 v);
+AliAnalysisTaskJetShapeTool(TClonesArray *list);
+ ~AliAnalysisTaskJetShapeTool();
+
+
+ void SetVecJ(TVector3 v);
+ void SetPtMinTr(Double_t a, Double_t b) {fPtMinTr[0] = a; fPtMinTr[1] = b;}
+ void SetNEntries(Int_t n) {fEntries = n;}
+
+ TArrayI GetListJ(Int_t b, Int_t i)  {return fListJ[b][i];}
+ TArrayI GetListB1(Int_t b, Int_t i) {return fListB1[b][i];}
+ TArrayI GetListB2(Int_t b, Int_t i) {return fListB2[b][i];}
+
+ Int_t GetSizeJ(Int_t b, Int_t i)  {return fListJc[b][i];}
+ Int_t GetSizeB1(Int_t b, Int_t i) {return fListB1c[b][i];}
+ Int_t GetSizeB2(Int_t b, Int_t i) {return fListB2c[b][i];}
+
+ TVector3 *GetJ(Int_t b,Int_t i, Int_t p)
+  { return (TVector3*)GetAt(fListJ[b][i].At(p));  }
+ TVector3 *GetB1(Int_t b,Int_t i, Int_t p)
+  { return (TVector3*)GetAt(fListB1[b][i].At(p)); }
+ TVector3 *GetB2(Int_t b,Int_t i, Int_t p)
+  { return (TVector3*)GetAt(fListB2[b][i].At(p)); }
+
+ TVector3 GetVecJ() {return  fvecJ;}
+ TVector3 GetVecB1() {return fvecB1;}
+ TVector3 GetVecB2() {return fvecB2;}
+
+ Double_t GetR(Int_t b) {return fR[b];}
+
+ TVector3 *GetAt(Int_t i)
+ { 
+   if(i<0 || i>= fEntries) printf(" TVector3 *GetAt(Int_t i) for i= %d\n",i);
+return (TVector3*) fList->At(i); 
+}
+ Double_t GetLocPhi(TVector3 v, Int_t i)
+ {
+   TVector3 p(*GetAt(i));
+   //   GetLocalMom(v, &vtmp);
+   //   return CalcdPhi0To2pi(vtmp.Phi());
+
+   Double_t phi =  TMath::ATan2(p.Eta() - v.Eta(), CalcdPhi(p.Phi(), v.Phi()) );
+
+   return CalcdPhi0To2pi(phi);
+ }
+
+
+ Double_t GetLocPhiJ(Int_t b, Int_t i, Int_t index)
+ {
+   TVector3 v1(*GetAt(fListJ[b][i].At(index)));
+   Double_t phi =  TMath::ATan2(v1.Eta() - fvecJ.Eta(), CalcdPhi(v1.Phi(), fvecJ.Phi()) );
+   return CalcdPhi0To2pi(phi);
+ }
+
+
+ Double_t GetLocPhiB1(Int_t b, Int_t i, Int_t index) 
+ {
+   TVector3 v0(*GetAt(fListB1[b][i].At(index)));
+   Double_t cV = (fvecB1(0)*fvecJ(0) + fvecB1(1)*fvecJ(1))/fvecJ.Perp2();
+   Double_t sV = (fvecB1(1)*fvecJ(0) - fvecB1(0)*fvecJ(1))/fvecJ.Perp2();
+   TVector3 v1(v0(0)*cV+v0(1)*sV, v0(1)*cV-v0(0)*sV, v0(2));
+   //   GetLocalMom(fvecJ, &v1);
+   //   return CalcdPhi0To2pi(v1.Phi());
+
+   Double_t phi =  TMath::ATan2(v1.Eta() - fvecJ.Eta(), CalcdPhi(v1.Phi(), fvecJ.Phi()) );
+   return CalcdPhi0To2pi(phi);
+
+ }
+ Double_t GetLocPhiB2(Int_t b, Int_t i, Int_t index) 
+ {
+   TVector3 v0(*GetAt(fListB2[b][i].At(index)));
+   Double_t cV = (fvecB2(0)*fvecJ(0) + fvecB2(1)*fvecJ(1))/fvecJ.Perp2();
+   Double_t sV = (fvecB2(1)*fvecJ(0) - fvecB2(0)*fvecJ(1))/fvecJ.Perp2();
+   TVector3 v1(v0(0)*cV+v0(1)*sV, v0(1)*cV-v0(0)*sV, v0(2));
+   //   GetLocalMom(fvecJ, &v1);
+   //   return CalcdPhi0To2pi(v1.Phi());
+   Double_t phi =  TMath::ATan2(v1.Eta() - fvecJ.Eta(), CalcdPhi(v1.Phi(), fvecJ.Phi()) );
+   return CalcdPhi0To2pi(phi);
+
+ }
+
+
+ void Add(TVector3 v) 
+   // { new(fList[fList.GetEntriesFast()]) TVector3(v);}
+ { new((*fList)[fEntries]) TVector3(v); fEntries++;}
+
+ void AddToJ(Int_t b, Int_t i, Int_t index) {
+   fListJ[b][i].AddAt(index, fListJc[b][i]);
+   fListJc[b][i]++;
+ }
+
+ void AddToB1(Int_t b, Int_t i, Int_t index) {
+   fListB1[b][i].AddAt(index, fListB1c[b][i]);
+   fListB1c[b][i]++;
+ }
+
+ void AddToB2(Int_t b, Int_t i, Int_t index) {
+   fListB2[b][i].AddAt(index, fListB2c[b][i]);
+   fListB2c[b][i]++;
+ }
+
+
+
+ TVector3 GetPRecJ(Int_t b, Int_t i) {return fPRecJ[b][i];}
+ TVector3 GetPRecInRJ() {return fPRecInRJ;}//
+
+
+ void Init();
+ // Bool_t FindCorrelationAxesAnd(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario=0);
+ Bool_t FindCorrelationAxes(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario=3);
+ Bool_t FindCorrelationAxesCorr(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario, Double_t R);
+
+ // virtual  void  Print(Option_t* /*option = ""*/) const;
+ void  PRINT() const;
+ void PRINT(TArrayI a, Int_t n, const char *txt="") const;
+
+ void Clean();
+
+ Double_t CalcdPhi0To2pi(Double_t phi1, Double_t phi2=0.)
+{
+  Double_t dphi = CalcdPhi(phi1, phi2);
+  while(dphi<0) dphi+=TMath::TwoPi();
+  while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
+  return dphi;
+}
+
+
+private:
+
+ TVector3 fvecJ;
+ TVector3 fvecB1;
+ TVector3 fvecB2;
+
+
+ Double_t fRmax;
+
+ static const Int_t fgkbinR= 3; // n bins
+ Double_t fR[fgkbinR]; 
+ TArrayI fListJ[fgkbinR][2]; //
+ TArrayI fListB1[fgkbinR][2];//
+ TArrayI fListB2[fgkbinR][2];//
+
+ Int_t fListJc[fgkbinR][2]; //
+ Int_t fListB1c[fgkbinR][2];//
+ Int_t fListB2c[fgkbinR][2];//
+
+ TVector3 fPRecJ[fgkbinR][2];//
+ TVector3 fPRecInRJ;//
+
+ Double_t fPtMinTr[2];
+
+
+ TClonesArray *fList; //!
+
+ Int_t fEntries;
+
+ Double_t CalcR(TVector3 v1, TVector3 v2);
+ Double_t CalcdPhi(Double_t phi1, Double_t phi2);
+ void GetLocalMom(TVector3 vecJ, TVector3 *q);
+
+
+  AliAnalysisTaskJetShapeTool(const AliAnalysisTaskJetShapeTool&);            // not implemented
+  AliAnalysisTaskJetShapeTool& operator=(const AliAnalysisTaskJetShapeTool&); // not implemented
+
+ ClassDef(AliAnalysisTaskJetShapeTool, 1)   // tbd
+
+};
+
+
+
+
+class AliAnalysisTaskJetShapeHM :public TObject {
+
+public:
+  AliAnalysisTaskJetShapeHM(TString comment = "", Bool_t kMC = kFALSE);
+ ~AliAnalysisTaskJetShapeHM();
+
+
+ void SetPtJetBins(Int_t Nbin=-1, Double_t *array = 0);
+ void SetRBins(    Int_t Nbin = 8, Double_t Rmax = 0.4) {fPsiVsRNbin = Nbin; fRmax = Rmax;}
+ void SetPhiNbins(Int_t Nbin = 1024) {fPsiVsPhiNbin = Nbin;}
+ void SetFilterMask(UInt_t i = 0){fFilterMask = i;}
+ void SetEtaTrackMax(Double_t e =0.9) {fEtaTrackMax = e;}
+ void SetPtTrackRange(Double_t min = 1., Double_t max = 100.) {fPtTrackMin = min; fPtTrackMax = max;}
+ void SetPtJetRange(Double_t min = 1., Double_t max = 200. ) {fPtJmin  =min; fPtJmax = max; }
+
+
+ Double_t GetPtJ() {return fPtJ;}
+ void InitHistos();
+ Bool_t AddEvent(AliAODEvent* aodE,  AliAODJet *jet, Double_t DeltaPtJ=0.);
+ void AddToList(TList *list);
+
+
+private:
+ TString fComment;
+ Bool_t fkMC;
+ TH1F *fhEvent;//! 
+ TH1F *fhMult; //! 
+ TH1F *fhPtJ;//! 
+ TH2F *fhPtJvsPtCorr;//! 
+ TH1F *fhPsiVsR;//! 
+ TH2F *fhPsiVsRPtJ;//! 
+ TH2F *fhPhiEtaTrack;//! 
+
+ TH1F *fhTMA_JAA[3];//! 
+ TH1F *fhTMA_JAp[3];//! 
+ TH1F *fhTMA_B1AA[3];//! 
+ TH1F *fhTMA_B1Ap[3];//! 
+ TH1F *fhTMA_B2AA[3];//! 
+ TH1F *fhTMA_B2Ap[3];//! 
+
+
+ Int_t     fPtJetNbin;
+ TArrayD   fPtJetArray;
+ Int_t     fPsiVsRNbin;
+ Double_t  fRmax;
+ Int_t     fPsiVsPhiNbin;
+
+ UInt_t   fFilterMask;
+ Double_t fEtaTrackMax;
+ Double_t fPtTrackMin;
+ Double_t fPtTrackMax;
+ Double_t fPtJmin;
+ Double_t fPtJmax;
+ Double_t fPtJ;
+
+ TVector3 fJvec;
+
+
+ Double_t CalcR(TVector3 v1, TVector3 v2);
+ Double_t CalcdPhi(Double_t phi1, Double_t phi2);
+
+ TH1F* Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax,  const char* xLabel = "", Int_t color=1, const char* yLabel="");
+ TH1F* Hist1D(const char* name, Int_t nBins, Double_t *xArray,  const char* xLabel = "", Int_t color=1, const char* yLabel="");
+ TH2F *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 = NULL, const char* yLabel = NULL, Int_t color=1);
+ TH2F *Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t *yArray, const char* xLabel = NULL, const char* yLabel = NULL, Int_t color=1, const char* zLabel = NULL);
+ TH2F *Hist2D(const char* name, Int_t nBinsx, Double_t *yArrax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel = NULL, const char* yLabel = NULL, Int_t color=1, const char* zLabel = NULL);
+
+
+  AliAnalysisTaskJetShapeHM(const AliAnalysisTaskJetShapeHM&);            // not implemented
+  AliAnalysisTaskJetShapeHM& operator=(const AliAnalysisTaskJetShapeHM&); // not implemented
+
+
+ClassDef(AliAnalysisTaskJetShapeHM,1)   // tbd
+};
+
+
+
+
+  AliAnalysisTaskJetShape(const char *name = "");
+  ~AliAnalysisTaskJetShape();
+
+  virtual void SetIsMC(Bool_t ismc=kTRUE) {fkMC=ismc;} 
+  virtual void SetCMSE(Double_t s = 7000.) {fCMSE=s;}
+
+
+  virtual void     SetNonStdFile(TString c){fNonStdFile = c;} 
+  virtual void     SetBranchNames(const TString &branch1, const TString &branch2);
+  virtual void     SetBackgroundBranch(TString &branch1, TString &branch2) { fBackgroundBranch[0] = branch1;  fBackgroundBranch[1] = branch2;};
+
+  virtual void     SetFilterMask(UInt_t i){fFilterMask = i;}
+  virtual void     SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
+  virtual void     SetCentMin(Float_t cent) { fCentMin = cent; }
+  virtual void     SetCentMax(Float_t cent) { fCentMax = cent; }
+  virtual void     SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; }
+  virtual void     SetEvtClassMax(Int_t evtClass) { fEvtClassMax = evtClass; }
+  virtual void     SetJetPtCorrMin(Float_t ptJ=20, Float_t ptB=20) { fPtJcorrMin = ptJ; fPtJBcorrMin = ptB;}
+
+  virtual void     SetPbPb(Bool_t a = kTRUE) {fkIsPbPb = a;}
+  // vtx
+  virtual void     SetVtxZRange(Double_t zmin=-10., Double_t zmax=10.) {fVtxZMin = zmin; fVtxZMax = zmax;}
+  virtual void     SetVixMinContrib(Int_t n=1) { fVtxMinContrib = n; }
+
+   //  virtual void     SetAngStructCloseTracks(Int_t yesno){fAngStructCloseTracks=yesno;}
+
+
+  //  virtual void   ConnectInputData(Option_t *);
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+
+
+
+
+ private:
+  //  TTree *fOutputTree; //!Output tree
+
+  TList       *fOutputList; // Output list
+  AliESDEvent *fESD;    // ESD object
+  AliAODEvent *fAODIn;    // AOD event
+  AliAODEvent *fAODOut;    // AOD event
+  AliAODExtension  *fAODExtension; //! where we take the jets from can be input or output AOD
+
+  //  AliTriggerAnalysis * fTriggerAnalysis; // trigger analysis object, to get the offline triggers
+
+
+   Bool_t   fkMC;         //is MC
+   Double_t fCMSE;      //cms energy
+   UInt_t fRunNb;       //run number
+   Bool_t fkIsPhysSel;  //tbd
+   TString       fNonStdFile; // delta AOD file
+   UInt_t  fFilterMask;            // filter bit for slecected tracks
+   AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline triggers to accept
+   Float_t fCentMin;     // lower bound on centrality
+   Float_t fCentMax;     // upper bound on centrality
+   Int_t   fEvtClassMin;         // lower bound on event class
+   Int_t   fEvtClassMax;         // upper bound on event class
+   Double_t fPtJcorrMin ;
+   Double_t fPtJBcorrMin;
+   Double_t fJpPtmin;
+   Double_t fJaPtmin;
+   Int_t fVtxMinContrib;
+   Double_t fVtxZMin;
+   Double_t fVtxZMax;
+   Bool_t fkIsPbPb;
+
+
+   TString fBackgroundBranch[2];  //tbd
+   TString fJetBranchName[2]; //  name of jet branches to compare
+
+
+   static const Int_t fgkbinNCent=7; // tbd
+  Double_t fgkbinCent[fgkbinNCent+1] ;// {0, 5, 10,  20, 40, 60, 80, 100}; 
+
+  TH1F *fhPtJL ;//!
+  TH1F *fhAreaJL ;//!
+
+  //  TClonesArray farray;
+    AliAnalysisTaskJetShapeHM *fanalJetSubStrHM;//!
+    AliAnalysisTaskJetShapeHM *fanalJetSubStrMCHM;//!
+
+   Bool_t IsGoodEvent();
+
+
+  Double_t CalcdPhi(Double_t phi1, Double_t phi2);
+  Double_t CalcdPhi0To2pi(Double_t phi1, Double_t phi2);
+
+  TH1F *Hist1D(const char* name, Int_t nBins , Double_t xMin, Double_t xMax, const char* xLabel = NULL, Int_t color=1);
+  TH2F *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 = NULL, const char* yLabel = NULL, Int_t color=1);
+  TH3F *Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t yMin, Double_t yMax, Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel = NULL, const char* yLabel = NULL, const char* zLabel = NULL, Int_t color=1);
+
+
+  //  AliAnalysisTaskJetShape(const AnalysisJetMain&); // not implemented
+  //  AliAnalysisTaskJetShape& operator=(const AnalysisJetMain&); // not implemented
+
+
+
+  AliAnalysisTaskJetShape(const AliAnalysisTaskJetShape&);            // not implemented
+  AliAnalysisTaskJetShape& operator=(const AliAnalysisTaskJetShape&); // not implemented
+
+
+  ClassDef(AliAnalysisTaskJetShape, 1)
+};
+
+
+
+
+
+
+
+
+
+#endif