--- /dev/null
+// 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&¢rality>=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
+
+
+
+
+
+
+
+
--- /dev/null
+#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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff