[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalQGTagging.cxx

//
// Jet QG tagging analysis task.
//
// Author: D. Caffarri, L. Cunqueiro 

#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THnSparse.h>
#include <TTree.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TProfile.h>
#include <TChain.h>
#include <TSystem.h>
#include <TFile.h>
#include <TKey.h>
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"
#include "TVector3.h"
#include "TVector2.h"

#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliEmcalParticle.h"
#include "AliMCEvent.h"
#include "AliGenPythiaEventHeader.h"
#include "AliAODMCHeader.h"
#include "AliMCEvent.h"
#include "AliAnalysisManager.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"
#include "AliStackPartonInfo.h"


#include "AliAODEvent.h"

#include "AliAnalysisTaskEmcalQGTagging.h"

using std::cout;
using std::endl;

ClassImp(AliAnalysisTaskEmcalQGTagging)

//________________________________________________________________________
AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalQGTagging", kTRUE),
  fContainer(0),
  fMinFractionShared(0),
  fJetShapeType(kData),
  fJetShapeSub(kNoSub),
  fShapesVar(0),
  fPtThreshold(-9999.),
  fRMatching(0.3),
  fPhiJetCorr6(0x0), 
  fPhiJetCorr7(0x0),
  fEtaJetCorr6(0x0),
  fEtaJetCorr7(0x0),
  fPtJetCorr(0x0),
  fPtJet(0x0),
  fTreeObservableTagging(0)

{
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fContainer(0),
  fMinFractionShared(0),
  fJetShapeType(kData),
  fJetShapeSub(kNoSub),
  fShapesVar(0),
  fPtThreshold(-9999.),
  fRMatching(0.3),
  fPhiJetCorr6(0x0), 
  fPhiJetCorr7(0x0),
  fEtaJetCorr6(0x0),
  fEtaJetCorr7(0x0),
  fPtJetCorr(0x0),
  fPtJet(0x0),
  fTreeObservableTagging(0)
  
{
  // Standard constructor.
  
  SetMakeGeneralHistograms(kTRUE);

  DefineOutput(1, TTree::Class());

}

//________________________________________________________________________
AliAnalysisTaskEmcalQGTagging::~AliAnalysisTaskEmcalQGTagging()
{
  // Destructor.
}

//________________________________________________________________________
 void AliAnalysisTaskEmcalQGTagging::UserCreateOutputObjects()
{
  // Create user output.

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
  Int_t nVar = 17; 
  fShapesVar = new Float_t [nVar]; 
  TString *fShapesVarNames = new TString [nVar];

  fShapesVarNames[0] = "partonCode"; 
  fShapesVarNames[1] = "ptJet"; 
  fShapesVarNames[2] = "ptDJet"; 
  fShapesVarNames[3] = "mJet";
  fShapesVarNames[4] = "nbOfConst";
  fShapesVarNames[5] = "angularity";
  fShapesVarNames[6] = "circularity";
  fShapesVarNames[7] = "lesub";
  fShapesVarNames[8] = "sigma2";

  fShapesVarNames[9] = "ptJetMatch"; 
  fShapesVarNames[10] = "ptDJetMatch"; 
  fShapesVarNames[11] = "mJetMatch";
  fShapesVarNames[12] = "nbOfConstMatch";
  fShapesVarNames[13] = "angularityMatch";
  fShapesVarNames[14] = "circularityMatch";
  fShapesVarNames[15] = "lesubMatch";
  fShapesVarNames[16] = "sigma2Match";
  
  for(Int_t ivar=0; ivar < nVar; ivar++){
    cout<<"looping over variables"<<endl;
    fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));

    //if( ivar == 4 )  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/I", fShapesVarNames[ivar].Data()));

  }
  
  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
  fOutput->Add(fPhiJetCorr6);
  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
  fOutput->Add(fEtaJetCorr6);
  
  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
  fOutput->Add(fPhiJetCorr7);
  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
  fOutput->Add(fEtaJetCorr7);
  
  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200,  100, 0, 200);
  fOutput->Add(fPtJetCorr);
  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
  fOutput->Add(fPtJet); 


  fOutput->Add(fTreeObservableTagging);
  TH1::AddDirectory(oldStatus);
  PostData(1, fOutput); // Post data for ALL output slots > 0 here.

}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalQGTagging::Run()
{
  // Run analysis code here, if needed. It will be executed before FillHistograms().

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalQGTagging::FillHistograms()
{
  // Fill histograms.
  //cout<<"base container"<<endl;
  AliEmcalJet* jet1 = NULL;
  AliJetContainer *jetCont = GetJetContainer(0);

  if(jetCont) {
    jetCont->ResetCurrentID();
    while((jet1 = jetCont->GetNextAcceptJet())) {
      if (!jet1) continue;
      AliEmcalJet* jet2 = 0x0;
      fPtJet->Fill(jet1->Pt());

      if (!(fJetShapeType == kData)) {
	AliStackPartonInfo *partonsInfo = 0x0;
	if((fJetShapeType == kTrueDet) || (fJetShapeType == kDetEmb)){
	  AliJetContainer *jetContTrue = GetJetContainer(1);
	  jet2 = jet1->ClosestJet();
	  if (!jet2) {
	    Printf("jet2 not exists, returning");
	    continue;
	  }
	  
	  Double_t fraction = jetCont->GetFractionSharedPt(jet1);
          cout<<"hey a jet"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<endl;
	  if(fraction<fMinFractionShared) continue;
	  partonsInfo = (AliStackPartonInfo*) jetContTrue->GetStackPartonInfo();     
	  if(!partonsInfo) return 0;
	}
	else {
	  partonsInfo = (AliStackPartonInfo*) jetCont->GetStackPartonInfo(); 
	  jet2=jet1;
          if(!partonsInfo) return 0;
	}
	
	Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6()); 
	Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
     	
	if (jp1< -1*TMath::Pi()) jp1 = (-2*TMath::Pi())-jp1;
	else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
	Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
	fEtaJetCorr6->Fill(jet2->Eta(), partonsInfo->GetPartonEta6());
	fPhiJetCorr6->Fill(jet2->Phi(), partonsInfo->GetPartonPhi6());
	if(dRp1 < fRMatching) {
	  fShapesVar[0] = partonsInfo->GetPartonFlag6();
	  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet2->Pt());
	}
	else {
	  jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi7());
	  detap1=(jet2->Eta())-(partonsInfo->GetPartonEta7());
	  if (jp1< -1*TMath::Pi()) jp1= (-2*TMath::Pi())-jp1;
	  else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
	  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
	  fEtaJetCorr7->Fill(jet2->Eta(), partonsInfo->GetPartonEta7());
	  fPhiJetCorr7->Fill(jet2->Phi(), partonsInfo->GetPartonPhi7());
	  if(dRp1 < fRMatching) {
	    fShapesVar[0] = partonsInfo->GetPartonFlag7();
	    fPtJetCorr ->Fill(partonsInfo->GetPartonPt7(), jet2->Pt());
	  }
	  else continue;
	}
      }
      else
	fShapesVar[0] = 0.;
     
      Double_t ptSubtracted = 0; 

      if (((fJetShapeType == kData) || (fJetShapeType == kDetEmb)) && (fJetShapeSub == kConstSub))
	  ptSubtracted = jet1->Pt(); 
      else ptSubtracted  = jet1->Pt() - GetRhoVal(0)*jet1->Area();
	  
      if ((fJetShapeType == kData || fJetShapeType== kDetEmb)) 
	if ( ptSubtracted < fPtThreshold) continue;	
      
      fShapesVar[1] = ptSubtracted;
      fShapesVar[2] = GetJetpTD(jet1,0);
      fShapesVar[3] = GetJetMass(jet1);
      fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1);
      fShapesVar[5] = GetJetAngularity(jet1,0);
      fShapesVar[6] = GetJetCircularity(jet1,0);
      fShapesVar[7] = GetJetLeSub(jet1,0);
      fShapesVar[8] = GetSigma2(jet1,0);
      
      Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0.;
      Int_t kMatched = 0;
      if (fJetShapeType == kTrueDet || fJetShapeType == kDetEmb) {
	kMatched = 1;
	ptMatch=jet2->Pt();
	ptDMatch=GetJetpTD(jet2, kMatched); 
	massMatch=GetJetMass(jet2);
	constMatch=1.*GetJetNumberOfConstituents(jet2);
	angulMatch=GetJetAngularity(jet2, kMatched);
	circMatch=GetJetCircularity(jet2, kMatched);
	lesubMatch=GetJetLeSub(jet2, kMatched);
	sigma2Match = GetSigma2(jet2, kMatched);
      }

      if (fJetShapeType == kTrue || fJetShapeType == kData) {
	kMatched = 0;
	ptMatch=0.;
	ptDMatch=0.; 
	massMatch=0.;
	constMatch=0.;
	angulMatch=0.;
	circMatch=0.;
	lesubMatch=0.;
	sigma2Match =0.;
      }
      
      fShapesVar[9] = ptMatch;
      fShapesVar[10] = ptDMatch;
      fShapesVar[11] = massMatch;
      fShapesVar[12] = constMatch;
      fShapesVar[13] = angulMatch;
      fShapesVar[14] = circMatch;
      fShapesVar[15] = lesubMatch;
      fShapesVar[16] = sigma2Match;

      fTreeObservableTagging->Fill();
       
    }
    
  } 
  
  return kTRUE;
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet) {
  //calc subtracted jet mass
  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtracted();
  else 
    return jet->M();
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){

  AliJetContainer *jetCont = GetJetContainer(jetContNb);
  if (!jet->GetNumberOfTracks())
      return 0; 
    Double_t den=0.;
    Double_t num = 0.;
    AliVParticle *vp1 = 0x0;
    for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
      vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
      Double_t dphi = vp1->Phi()-jet->Phi();
      if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
      if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
      Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
      Double_t dr = TMath::Sqrt(dr2);
      num=num+vp1->Pt()*dr;
      den=den+vp1->Pt();
    }
    return num/den;
} 

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) {

  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtractedAngularity();
  else
    return Angularity(jet, jetContNb);
 
}


//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){

  AliJetContainer *jetCont = GetJetContainer(jetContNb);
  if (!jet->GetNumberOfTracks())
      return 0; 
    Double_t den=0.;
    Double_t num = 0.;
    AliVParticle *vp1 = 0x0;
    for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
      vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
      num=num+vp1->Pt()*vp1->Pt();
      den=den+vp1->Pt();
    }
    return TMath::Sqrt(num)/den;
} 

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) {
  //calc subtracted jet mass
  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtractedpTD();
  else
    return PTD(jet, jetContNb);
 
}

//_____________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){

  AliJetContainer *jetCont = GetJetContainer(jetContNb);
  if (!jet->GetNumberOfTracks())
    return 0; 
  Double_t mxx    = 0.;
  Double_t myy    = 0.;
  Double_t mxy    = 0.;
  int  nc     = 0;
  Double_t sump2  = 0.;
  Double_t pxjet=jet->Px();
  Double_t pyjet=jet->Py();
  Double_t pzjet=jet->Pz();
  
  
  //2 general normalized vectors perpendicular to the jet
  TVector3  ppJ1(pxjet, pyjet, pzjet);
  TVector3  ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
  ppJ3.SetMag(1.);
  TVector3  ppJ2(-pyjet, pxjet, 0);
  ppJ2.SetMag(1.);
  AliVParticle *vp1 = 0x0;
  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
    vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
    
    
    TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
   
    //local frame
    TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
    TVector3 pPerp = pp - pLong;
    //projection onto the two perpendicular vectors defined above
    
    Float_t ppjX = pPerp.Dot(ppJ2);
    Float_t ppjY = pPerp.Dot(ppJ3);
    Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
    if(ppjT<=0) return 0;
    
    mxx += (ppjX * ppjX / ppjT);
    myy += (ppjY * ppjY / ppjT);
    mxy += (ppjX * ppjY / ppjT);
    nc++;
    sump2 += ppjT;}
  
  if(nc<2) return 0;
  if(sump2==0) return 0;
  // Sphericity Matrix
  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 jev = 0;
  //  cout<<eval[0]<<" "<<eval[1]<<endl;
  if (eval[0] < eval[1]) jev = 1;
  TVectorD evec0(2);
  // Principle axis
  evec0 = TMatrixDColumn(evecm, jev);
  Double_t compx=evec0[0];
  Double_t compy=evec0[1];
  TVector2 evec(compx, compy);
  Double_t circ=0;
  if(jev==1) circ=2*eval[0];
  if(jev==0) circ=2*eval[1];
  
  return circ;
  
  
}


//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) {
  //calc subtracted jet mass
 
  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtractedCircularity();
  else
    return Circularity(jet, jetContNb);
 
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){

  AliJetContainer *jetCont = GetJetContainer(jetContNb);
  if (!jet->GetNumberOfTracks())
      return 0; 
    Double_t den=0.;
    Double_t num = 0.;
    AliVParticle *vp1 = 0x0;
    AliVParticle *vp2 = 0x0;
    std::vector<int> ordindex;
    ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
    
   vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));  
   vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));  
     
  num=vp1->Pt();
  den=vp2->Pt();
  
return num-den;
} 

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
  //calc subtracted jet mass
 
  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtractedLeSub();
  else
    return LeSub(jet, jetContNb);
 
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet) {
  //calc subtracted jet mass
  
  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtractedConstituent();
  else
    return jet->GetNumberOfTracks();
 
}
   

//______________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){

  AliJetContainer *jetCont = GetJetContainer(jetContNb);
  if (!jet->GetNumberOfTracks())
      return 0; 
      Double_t mxx    = 0.;
      Double_t myy    = 0.;
      Double_t mxy    = 0.;
      int  nc     = 0;
      Double_t sump2  = 0.;
       
     AliVParticle *vp1 = 0x0;
     for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
       vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
       Double_t ppt=vp1->Pt();
       Double_t dphi = vp1->Phi()-jet->Phi();
       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
       Double_t deta = vp1->Eta()-jet->Eta();
       mxx += ppt*ppt*deta*deta;
       myy += ppt*ppt*dphi*dphi;
       mxy -= ppt*ppt*deta*dphi;
       nc++;
       sump2 += ppt*ppt;
       
     }  
     if(nc<2) return 0;
     if(sump2==0) return 0;
     // Sphericity Matrix
     Double_t ele[4] = {mxx , mxy , mxy , myy };
     TMatrixDSym m0(2,ele);
     
     // Find eigenvectors
     TMatrixDSymEigen m(m0);
     TVectorD eval(2);
     TMatrixD evecm = m.GetEigenVectors();
     eval  = m.GetEigenValues();
     // Largest eigenvector
     int jev = 0;
     //  cout<<eval[0]<<" "<<eval[1]<<endl;
     if (eval[0] < eval[1]) jev = 1;
     TVectorD evec0(2);
     // Principle axis
     evec0 = TMatrixDColumn(evecm, jev);
     Double_t compx=evec0[0];
     Double_t compy=evec0[1];
     TVector2 evec(compx, compy);
     Double_t sig=0;
     if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
     if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
     
     return sig;
     
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) {
  //calc subtracted jet mass
 
  if(fJetShapeSub==kDerivSub)
    return jet->GetSecondOrderSubtractedSigma2();
  else
    return Sigma2(jet, jetContNb);
 
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
  //
  // retrieve event objects
  //
  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;

  return kTRUE;
}

//_______________________________________________________________________
void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.

  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
  // if (!fTreeObservableTagging){
  //   Printf("ERROR: fTreeObservableTagging not available"); 
  //   return;
  // }

}
Commit	Line	Data
cfa7b4df	1	//
	2	// Jet QG tagging analysis task.
	3	//
	4	// Author: D. Caffarri, L. Cunqueiro
	5
	6	#include <TClonesArray.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TH3F.h>
	10	#include <THnSparse.h>
	11	#include <TTree.h>
	12	#include <TList.h>
	13	#include <TLorentzVector.h>
	14	#include <TProfile.h>
	15	#include <TChain.h>
	16	#include <TSystem.h>
	17	#include <TFile.h>
	18	#include <TKey.h>
	19	#include "TMatrixD.h"
	20	#include "TMatrixDSym.h"
	21	#include "TMatrixDSymEigen.h"
	22	#include "TVector3.h"
	23	#include "TVector2.h"
	24
	25	#include "AliVCluster.h"
	26	#include "AliVTrack.h"
	27	#include "AliEmcalJet.h"
	28	#include "AliRhoParameter.h"
	29	#include "AliLog.h"
	30	#include "AliEmcalParticle.h"
	31	#include "AliMCEvent.h"
	32	#include "AliGenPythiaEventHeader.h"
	33	#include "AliAODMCHeader.h"
	34	#include "AliMCEvent.h"
	35	#include "AliAnalysisManager.h"
	36	#include "AliJetContainer.h"
	37	#include "AliParticleContainer.h"
	38	#include "AliStackPartonInfo.h"
	39
	40
	41	#include "AliAODEvent.h"
	42
	43	#include "AliAnalysisTaskEmcalQGTagging.h"
	44
e8c28d9e	45	using std::cout;
	46	using std::endl;
	47
cfa7b4df	48	ClassImp(AliAnalysisTaskEmcalQGTagging)
	49
	50	//________________________________________________________________________
	51	AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging() :
	52	AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalQGTagging", kTRUE),
	53	fContainer(0),
	54	fMinFractionShared(0),
a68b1989	55	fJetShapeType(kData),
a68b1989	56	fJetShapeSub(kNoSub),
cfa7b4df	57	fShapesVar(0),
34b1e79e	58	fPtThreshold(-9999.),
cfa7b4df	59	fRMatching(0.3),
	60	fPhiJetCorr6(0x0),
	61	fPhiJetCorr7(0x0),
	62	fEtaJetCorr6(0x0),
	63	fEtaJetCorr7(0x0),
	64	fPtJetCorr(0x0),
	65	fPtJet(0x0),
	66	fTreeObservableTagging(0)
34b1e79e	67
cfa7b4df	68	{
	69	SetMakeGeneralHistograms(kTRUE);
	70	}
	71
	72	//________________________________________________________________________
	73	AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) :
	74	AliAnalysisTaskEmcalJet(name, kTRUE),
	75	fContainer(0),
	76	fMinFractionShared(0),
a68b1989	77	fJetShapeType(kData),
a68b1989	78	fJetShapeSub(kNoSub),
cfa7b4df	79	fShapesVar(0),
34b1e79e	80	fPtThreshold(-9999.),
cfa7b4df	81	fRMatching(0.3),
	82	fPhiJetCorr6(0x0),
	83	fPhiJetCorr7(0x0),
	84	fEtaJetCorr6(0x0),
	85	fEtaJetCorr7(0x0),
	86	fPtJetCorr(0x0),
	87	fPtJet(0x0),
	88	fTreeObservableTagging(0)
34b1e79e	89
cfa7b4df	90	{
	91	// Standard constructor.
	92
	93	SetMakeGeneralHistograms(kTRUE);
	94
	95	DefineOutput(1, TTree::Class());
	96
	97	}
	98
	99	//________________________________________________________________________
	100	AliAnalysisTaskEmcalQGTagging::~AliAnalysisTaskEmcalQGTagging()
	101	{
	102	// Destructor.
	103	}
	104
	105	//________________________________________________________________________
	106	void AliAnalysisTaskEmcalQGTagging::UserCreateOutputObjects()
	107	{
	108	// Create user output.
	109
	110	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
	111
	112	Bool_t oldStatus = TH1::AddDirectoryStatus();
	113	TH1::AddDirectory(kFALSE);
	114
	115	fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
a68b1989	116	Int_t nVar = 17;
cfa7b4df	117	fShapesVar = new Float_t [nVar];
	118	TString *fShapesVarNames = new TString [nVar];
	119
	120	fShapesVarNames[0] = "partonCode";
	121	fShapesVarNames[1] = "ptJet";
	122	fShapesVarNames[2] = "ptDJet";
	123	fShapesVarNames[3] = "mJet";
	124	fShapesVarNames[4] = "nbOfConst";
	125	fShapesVarNames[5] = "angularity";
	126	fShapesVarNames[6] = "circularity";
	127	fShapesVarNames[7] = "lesub";
	128	fShapesVarNames[8] = "sigma2";
a68b1989	129
	130	fShapesVarNames[9] = "ptJetMatch";
	131	fShapesVarNames[10] = "ptDJetMatch";
	132	fShapesVarNames[11] = "mJetMatch";
	133	fShapesVarNames[12] = "nbOfConstMatch";
	134	fShapesVarNames[13] = "angularityMatch";
	135	fShapesVarNames[14] = "circularityMatch";
	136	fShapesVarNames[15] = "lesubMatch";
	137	fShapesVarNames[16] = "sigma2Match";
cfa7b4df	138
	139	for(Int_t ivar=0; ivar < nVar; ivar++){
	140	cout<<"looping over variables"<<endl;
	141	fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));
	142
	143	//if( ivar == 4 ) fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/I", fShapesVarNames[ivar].Data()));
	144
	145	}
	146
	147	fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2TMath::Pi(), 50, 0, 2TMath::Pi());
	148	fOutput->Add(fPhiJetCorr6);
	149	fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
	150	fOutput->Add(fEtaJetCorr6);
	151
	152	fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2TMath::Pi(), 50, 0, 2TMath::Pi());
	153	fOutput->Add(fPhiJetCorr7);
	154	fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
	155	fOutput->Add(fEtaJetCorr7);
	156
	157	fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
	158	fOutput->Add(fPtJetCorr);
	159	fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
	160	fOutput->Add(fPtJet);
	161
	162
	163	fOutput->Add(fTreeObservableTagging);
	164	TH1::AddDirectory(oldStatus);
	165	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	166
	167	}
	168
	169	//________________________________________________________________________
	170	Bool_t AliAnalysisTaskEmcalQGTagging::Run()
	171	{
	172	// Run analysis code here, if needed. It will be executed before FillHistograms().
	173
	174	return kTRUE;
	175	}
	176
	177	//________________________________________________________________________
	178	Bool_t AliAnalysisTaskEmcalQGTagging::FillHistograms()
	179	{
	180	// Fill histograms.
	181	//cout<<"base container"<<endl;
	182	AliEmcalJet* jet1 = NULL;
	183	AliJetContainer *jetCont = GetJetContainer(0);
a68b1989	184
cfa7b4df	185	if(jetCont) {
	186	jetCont->ResetCurrentID();
	187	while((jet1 = jetCont->GetNextAcceptJet())) {
	188	if (!jet1) continue;
a68b1989	189	AliEmcalJet* jet2 = 0x0;
cfa7b4df	190	fPtJet->Fill(jet1->Pt());
a68b1989	191
a68b1989	192	if (!(fJetShapeType == kData)) {
cfa7b4df	193	AliStackPartonInfo *partonsInfo = 0x0;
a68b1989	194	if((fJetShapeType == kTrueDet) \|\| (fJetShapeType == kDetEmb)){
cfa7b4df	195	AliJetContainer *jetContTrue = GetJetContainer(1);
	196	jet2 = jet1->ClosestJet();
	197	if (!jet2) {
	198	Printf("jet2 not exists, returning");
	199	continue;
	200	}
a68b1989	201
cfa7b4df	202	Double_t fraction = jetCont->GetFractionSharedPt(jet1);
a68b1989	203	cout<<"hey a jet"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<endl;
cfa7b4df	204	if(fraction<fMinFractionShared) continue;
a68b1989	205	partonsInfo = (AliStackPartonInfo*) jetContTrue->GetStackPartonInfo();
a68b1989	206	if(!partonsInfo) return 0;
cfa7b4df	207	}
	208	else {
	209	partonsInfo = (AliStackPartonInfo*) jetCont->GetStackPartonInfo();
	210	jet2=jet1;
be58fab7	211	if(!partonsInfo) return 0;
cfa7b4df	212	}
	213
	214	Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6());
	215	Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
	216
	217	if (jp1< -1TMath::Pi()) jp1 = (-2TMath::Pi())-jp1;
	218	else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
	219	Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
	220	fEtaJetCorr6->Fill(jet2->Eta(), partonsInfo->GetPartonEta6());
	221	fPhiJetCorr6->Fill(jet2->Phi(), partonsInfo->GetPartonPhi6());
	222	if(dRp1 < fRMatching) {
	223	fShapesVar[0] = partonsInfo->GetPartonFlag6();
	224	fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet2->Pt());
	225	}
	226	else {
	227	jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi7());
	228	detap1=(jet2->Eta())-(partonsInfo->GetPartonEta7());
	229	if (jp1< -1TMath::Pi()) jp1= (-2TMath::Pi())-jp1;
	230	else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
	231	dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
	232	fEtaJetCorr7->Fill(jet2->Eta(), partonsInfo->GetPartonEta7());
	233	fPhiJetCorr7->Fill(jet2->Phi(), partonsInfo->GetPartonPhi7());
	234	if(dRp1 < fRMatching) {
	235	fShapesVar[0] = partonsInfo->GetPartonFlag7();
	236	fPtJetCorr ->Fill(partonsInfo->GetPartonPt7(), jet2->Pt());
	237	}
	238	else continue;
	239	}
	240	}
	241	else
	242	fShapesVar[0] = 0.;
	243
34b1e79e	244	Double_t ptSubtracted = 0;
cfa7b4df	245
a68b1989	246	if (((fJetShapeType == kData) \|\| (fJetShapeType == kDetEmb)) && (fJetShapeSub == kConstSub))
	247	ptSubtracted = jet1->Pt();
	248	else ptSubtracted = jet1->Pt() - GetRhoVal(0)*jet1->Area();
	249
	250	if ((fJetShapeType == kData \|\| fJetShapeType== kDetEmb))
	251	if ( ptSubtracted < fPtThreshold) continue;
	252
34b1e79e	253	fShapesVar[1] = ptSubtracted;
a68b1989	254	fShapesVar[2] = GetJetpTD(jet1,0);
cfa7b4df	255	fShapesVar[3] = GetJetMass(jet1);
cfa7b4df	256	fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1);
a68b1989	257	fShapesVar[5] = GetJetAngularity(jet1,0);
	258	fShapesVar[6] = GetJetCircularity(jet1,0);
	259	fShapesVar[7] = GetJetLeSub(jet1,0);
	260	fShapesVar[8] = GetSigma2(jet1,0);
	261
	262	Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0.;
	263	Int_t kMatched = 0;
	264	if (fJetShapeType == kTrueDet \|\| fJetShapeType == kDetEmb) {
	265	kMatched = 1;
	266	ptMatch=jet2->Pt();
	267	ptDMatch=GetJetpTD(jet2, kMatched);
	268	massMatch=GetJetMass(jet2);
	269	constMatch=1.*GetJetNumberOfConstituents(jet2);
	270	angulMatch=GetJetAngularity(jet2, kMatched);
	271	circMatch=GetJetCircularity(jet2, kMatched);
	272	lesubMatch=GetJetLeSub(jet2, kMatched);
	273	sigma2Match = GetSigma2(jet2, kMatched);
	274	}
	275
	276	if (fJetShapeType == kTrue \|\| fJetShapeType == kData) {
	277	kMatched = 0;
	278	ptMatch=0.;
	279	ptDMatch=0.;
	280	massMatch=0.;
	281	constMatch=0.;
	282	angulMatch=0.;
	283	circMatch=0.;
	284	lesubMatch=0.;
	285	sigma2Match =0.;
	286	}
	287
	288	fShapesVar[9] = ptMatch;
	289	fShapesVar[10] = ptDMatch;
	290	fShapesVar[11] = massMatch;
	291	fShapesVar[12] = constMatch;
	292	fShapesVar[13] = angulMatch;
	293	fShapesVar[14] = circMatch;
	294	fShapesVar[15] = lesubMatch;
	295	fShapesVar[16] = sigma2Match;
	296
cfa7b4df	297	fTreeObservableTagging->Fill();
	298
	299	}
	300
	301	}
	302
	303	return kTRUE;
	304	}
	305
	306	//________________________________________________________________________
	307	Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet) {
	308	//calc subtracted jet mass
a68b1989	309	if(fJetShapeSub==kDerivSub)
cfa7b4df	310	return jet->GetSecondOrderSubtracted();
	311	else
	312	return jet->M();
	313	}
	314
	315	//________________________________________________________________________
a68b1989	316	Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
cfa7b4df	317
a68b1989	318	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	319	if (!jet->GetNumberOfTracks())
	320	return 0;
	321	Double_t den=0.;
	322	Double_t num = 0.;
	323	AliVParticle *vp1 = 0x0;
	324	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	325	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	326	Double_t dphi = vp1->Phi()-jet->Phi();
	327	if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
	328	if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
	329	Double_t dr2 = (vp1->Eta()-jet->Eta())(vp1->Eta()-jet->Eta()) + dphidphi;
	330	Double_t dr = TMath::Sqrt(dr2);
	331	num=num+vp1->Pt()*dr;
	332	den=den+vp1->Pt();
	333	}
	334	return num/den;
	335	}
	336
	337	//________________________________________________________________________
a68b1989	338	Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) {
cfa7b4df	339
a68b1989	340	if(fJetShapeSub==kDerivSub)
cfa7b4df	341	return jet->GetSecondOrderSubtractedAngularity();
cfa7b4df	342	else
a68b1989	343	return Angularity(jet, jetContNb);
cfa7b4df	344
	345	}
	346
	347
	348	//________________________________________________________________________
a68b1989	349	Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
cfa7b4df	350
a68b1989	351	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	352	if (!jet->GetNumberOfTracks())
	353	return 0;
	354	Double_t den=0.;
	355	Double_t num = 0.;
	356	AliVParticle *vp1 = 0x0;
	357	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	358	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	359	num=num+vp1->Pt()*vp1->Pt();
	360	den=den+vp1->Pt();
	361	}
	362	return TMath::Sqrt(num)/den;
	363	}
	364
	365	//________________________________________________________________________
a68b1989	366	Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) {
cfa7b4df	367	//calc subtracted jet mass
a68b1989	368	if(fJetShapeSub==kDerivSub)
cfa7b4df	369	return jet->GetSecondOrderSubtractedpTD();
cfa7b4df	370	else
a68b1989	371	return PTD(jet, jetContNb);
cfa7b4df	372
	373	}
	374
	375	//_____________________________________________________________________________
a68b1989	376	Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
cfa7b4df	377
a68b1989	378	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	379	if (!jet->GetNumberOfTracks())
	380	return 0;
	381	Double_t mxx = 0.;
	382	Double_t myy = 0.;
	383	Double_t mxy = 0.;
	384	int nc = 0;
	385	Double_t sump2 = 0.;
	386	Double_t pxjet=jet->Px();
	387	Double_t pyjet=jet->Py();
	388	Double_t pzjet=jet->Pz();
	389
	390
	391	//2 general normalized vectors perpendicular to the jet
	392	TVector3 ppJ1(pxjet, pyjet, pzjet);
	393	TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
	394	ppJ3.SetMag(1.);
	395	TVector3 ppJ2(-pyjet, pxjet, 0);
	396	ppJ2.SetMag(1.);
	397	AliVParticle *vp1 = 0x0;
	398	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	399	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	400
	401
	402	TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
	403
	404	//local frame
	405	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	406	TVector3 pPerp = pp - pLong;
	407	//projection onto the two perpendicular vectors defined above
	408
	409	Float_t ppjX = pPerp.Dot(ppJ2);
	410	Float_t ppjY = pPerp.Dot(ppJ3);
	411	Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
	412	if(ppjT<=0) return 0;
	413
	414	mxx += (ppjX * ppjX / ppjT);
	415	myy += (ppjY * ppjY / ppjT);
	416	mxy += (ppjX * ppjY / ppjT);
	417	nc++;
	418	sump2 += ppjT;}
	419
	420	if(nc<2) return 0;
	421	if(sump2==0) return 0;
	422	// Sphericity Matrix
	423	Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
	424	TMatrixDSym m0(2,ele);
	425
	426	// Find eigenvectors
	427	TMatrixDSymEigen m(m0);
	428	TVectorD eval(2);
	429	TMatrixD evecm = m.GetEigenVectors();
	430	eval = m.GetEigenValues();
	431	// Largest eigenvector
	432	int jev = 0;
	433	// cout<<eval[0]<<" "<<eval[1]<<endl;
	434	if (eval[0] < eval[1]) jev = 1;
	435	TVectorD evec0(2);
	436	// Principle axis
	437	evec0 = TMatrixDColumn(evecm, jev);
	438	Double_t compx=evec0[0];
	439	Double_t compy=evec0[1];
	440	TVector2 evec(compx, compy);
	441	Double_t circ=0;
	442	if(jev==1) circ=2*eval[0];
443	if(jev==0) circ=2*eval[1];
444
445	return circ;
446
447
448
449	}
450
451
452
453
454	//________________________________________________________________________
a68b1989	455	Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) {
cfa7b4df	456	//calc subtracted jet mass
cfa7b4df	457
a68b1989	458	if(fJetShapeSub==kDerivSub)
cfa7b4df	459	return jet->GetSecondOrderSubtractedCircularity();
cfa7b4df	460	else
a68b1989	461	return Circularity(jet, jetContNb);
cfa7b4df	462
	463	}
	464
	465	//________________________________________________________________________
a68b1989	466	Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
cfa7b4df	467
a68b1989	468	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	469	if (!jet->GetNumberOfTracks())
	470	return 0;
	471	Double_t den=0.;
	472	Double_t num = 0.;
	473	AliVParticle *vp1 = 0x0;
	474	AliVParticle *vp2 = 0x0;
	475	std::vector<int> ordindex;
	476	ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
	477
	478	vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
	479	vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
	480
	481	num=vp1->Pt();
	482	den=vp2->Pt();
	483
	484	return num-den;
	485	}
	486
	487	//________________________________________________________________________
a68b1989	488	Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
cfa7b4df	489	//calc subtracted jet mass
cfa7b4df	490
a68b1989	491	if(fJetShapeSub==kDerivSub)
cfa7b4df	492	return jet->GetSecondOrderSubtractedLeSub();
cfa7b4df	493	else
a68b1989	494	return LeSub(jet, jetContNb);
cfa7b4df	495
	496	}
	497
	498	//________________________________________________________________________
	499	Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet) {
	500	//calc subtracted jet mass
	501
a68b1989	502	if(fJetShapeSub==kDerivSub)
cfa7b4df	503	return jet->GetSecondOrderSubtractedConstituent();
	504	else
	505	return jet->GetNumberOfTracks();
	506
	507	}
	508
	509
	510	//______________________________________________________________________________
a68b1989	511	Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
cfa7b4df	512
a68b1989	513	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	514	if (!jet->GetNumberOfTracks())
	515	return 0;
	516	Double_t mxx = 0.;
	517	Double_t myy = 0.;
	518	Double_t mxy = 0.;
	519	int nc = 0;
	520	Double_t sump2 = 0.;
	521
	522	AliVParticle *vp1 = 0x0;
	523	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	524	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	525	Double_t ppt=vp1->Pt();
	526	Double_t dphi = vp1->Phi()-jet->Phi();
	527	if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
	528	if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
	529	Double_t deta = vp1->Eta()-jet->Eta();
	530	mxx += pptpptdeta*deta;
	531	myy += pptpptdphi*dphi;
	532	mxy -= pptpptdeta*dphi;
	533	nc++;
	534	sump2 += ppt*ppt;
	535
	536	}
	537	if(nc<2) return 0;
	538	if(sump2==0) return 0;
	539	// Sphericity Matrix
	540	Double_t ele[4] = {mxx , mxy , mxy , myy };
	541	TMatrixDSym m0(2,ele);
	542
	543	// Find eigenvectors
	544	TMatrixDSymEigen m(m0);
	545	TVectorD eval(2);
	546	TMatrixD evecm = m.GetEigenVectors();
	547	eval = m.GetEigenValues();
	548	// Largest eigenvector
	549	int jev = 0;
	550	// cout<<eval[0]<<" "<<eval[1]<<endl;
	551	if (eval[0] < eval[1]) jev = 1;
	552	TVectorD evec0(2);
	553	// Principle axis
	554	evec0 = TMatrixDColumn(evecm, jev);
	555	Double_t compx=evec0[0];
	556	Double_t compy=evec0[1];
	557	TVector2 evec(compx, compy);
	558	Double_t sig=0;
	559	if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
	560	if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
	561
	562	return sig;
	563
	564	}
	565
	566	//________________________________________________________________________
a68b1989	567	Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) {
cfa7b4df	568	//calc subtracted jet mass
cfa7b4df	569
a68b1989	570	if(fJetShapeSub==kDerivSub)
cfa7b4df	571	return jet->GetSecondOrderSubtractedSigma2();
cfa7b4df	572	else
a68b1989	573	return Sigma2(jet, jetContNb);
cfa7b4df	574
	575	}
	576
	577	//________________________________________________________________________
	578	Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
	579	//
	580	// retrieve event objects
	581	//
	582	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	583	return kFALSE;
	584
	585	return kTRUE;
	586	}
	587
	588	//_______________________________________________________________________
	589	void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *)
	590	{
	591	// Called once at the end of the analysis.
	592
	593	// fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
	594	// if (!fTreeObservableTagging){
	595	// Printf("ERROR: fTreeObservableTagging not available");
	596	// return;
	597	// }
	598
	599	}
	600