[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 "AliPythiaInfo.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);
  Float_t kWeight=1;
  if(jetCont) {
    jetCont->ResetCurrentID();
    while((jet1 = jetCont->GetNextAcceptJet())) {
      if (!jet1) continue;
      AliEmcalJet* jet2 = 0x0;
      fPtJet->Fill(jet1->Pt());

      if (!(fJetShapeType == kData)) {
	AliPythiaInfo *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 = (AliPythiaInfo*) jetContTrue->GetPythiaInfo();     
	  if(!partonsInfo) return 0;
          
	}
	else {
	  partonsInfo = (AliPythiaInfo*) jetCont->GetPythiaInfo(); 
	  jet2=jet1;
          if(!partonsInfo) return 0;
	}
	
	Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6()); 
	Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
     	kWeight=partonsInfo->GetPythiaEventWeight();
	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,0);
      fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
      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,kMatched);
	constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
	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();
      fTreeObservableTagging->SetWeight(kWeight); 
    }
    
  } 
  
  return kTRUE;
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0) {
  //calc subtracted jet mass
  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
    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) && (jetContNb==0))
    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)&&(jetContNb==0))
    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)&&(jetContNb==0))
    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)&&(jetContNb==0))
    return jet->GetSecondOrderSubtractedLeSub();
  else
    return LeSub(jet, jetContNb);
 
}

//________________________________________________________________________
Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0) {
  //calc subtracted jet mass
  
  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
    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)&&(jetContNb==0))
    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"
b929303d	38	#include "AliPythiaInfo.h"
cfa7b4df	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);
b929303d	184	Float_t kWeight=1;
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)) {
b929303d	193	AliPythiaInfo *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;
b929303d	205	partonsInfo = (AliPythiaInfo*) jetContTrue->GetPythiaInfo();
a68b1989	206	if(!partonsInfo) return 0;
b929303d	207
cfa7b4df	208	}
cfa7b4df	209	else {
b929303d	210	partonsInfo = (AliPythiaInfo*) jetCont->GetPythiaInfo();
cfa7b4df	211	jet2=jet1;
be58fab7	212	if(!partonsInfo) return 0;
cfa7b4df	213	}
	214
	215	Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6());
	216	Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
b929303d	217	kWeight=partonsInfo->GetPythiaEventWeight();
cfa7b4df	218	if (jp1< -1TMath::Pi()) jp1 = (-2TMath::Pi())-jp1;
	219	else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
	220	Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
	221	fEtaJetCorr6->Fill(jet2->Eta(), partonsInfo->GetPartonEta6());
	222	fPhiJetCorr6->Fill(jet2->Phi(), partonsInfo->GetPartonPhi6());
	223	if(dRp1 < fRMatching) {
	224	fShapesVar[0] = partonsInfo->GetPartonFlag6();
	225	fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet2->Pt());
	226	}
	227	else {
	228	jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi7());
	229	detap1=(jet2->Eta())-(partonsInfo->GetPartonEta7());
	230	if (jp1< -1TMath::Pi()) jp1= (-2TMath::Pi())-jp1;
	231	else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
	232	dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
	233	fEtaJetCorr7->Fill(jet2->Eta(), partonsInfo->GetPartonEta7());
	234	fPhiJetCorr7->Fill(jet2->Phi(), partonsInfo->GetPartonPhi7());
	235	if(dRp1 < fRMatching) {
	236	fShapesVar[0] = partonsInfo->GetPartonFlag7();
	237	fPtJetCorr ->Fill(partonsInfo->GetPartonPt7(), jet2->Pt());
	238	}
	239	else continue;
	240	}
	241	}
	242	else
	243	fShapesVar[0] = 0.;
	244
34b1e79e	245	Double_t ptSubtracted = 0;
cfa7b4df	246
a68b1989	247	if (((fJetShapeType == kData) \|\| (fJetShapeType == kDetEmb)) && (fJetShapeSub == kConstSub))
	248	ptSubtracted = jet1->Pt();
	249	else ptSubtracted = jet1->Pt() - GetRhoVal(0)*jet1->Area();
	250
	251	if ((fJetShapeType == kData \|\| fJetShapeType== kDetEmb))
	252	if ( ptSubtracted < fPtThreshold) continue;
	253
34b1e79e	254	fShapesVar[1] = ptSubtracted;
a68b1989	255	fShapesVar[2] = GetJetpTD(jet1,0);
e5aca8ab	256	fShapesVar[3] = GetJetMass(jet1,0);
e5aca8ab	257	fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
a68b1989	258	fShapesVar[5] = GetJetAngularity(jet1,0);
	259	fShapesVar[6] = GetJetCircularity(jet1,0);
	260	fShapesVar[7] = GetJetLeSub(jet1,0);
	261	fShapesVar[8] = GetSigma2(jet1,0);
	262
	263	Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0.;
	264	Int_t kMatched = 0;
	265	if (fJetShapeType == kTrueDet \|\| fJetShapeType == kDetEmb) {
	266	kMatched = 1;
	267	ptMatch=jet2->Pt();
	268	ptDMatch=GetJetpTD(jet2, kMatched);
e5aca8ab	269	massMatch=GetJetMass(jet2,kMatched);
e5aca8ab	270	constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
a68b1989	271	angulMatch=GetJetAngularity(jet2, kMatched);
	272	circMatch=GetJetCircularity(jet2, kMatched);
	273	lesubMatch=GetJetLeSub(jet2, kMatched);
	274	sigma2Match = GetSigma2(jet2, kMatched);
	275	}
	276
	277	if (fJetShapeType == kTrue \|\| fJetShapeType == kData) {
	278	kMatched = 0;
	279	ptMatch=0.;
	280	ptDMatch=0.;
	281	massMatch=0.;
	282	constMatch=0.;
	283	angulMatch=0.;
	284	circMatch=0.;
	285	lesubMatch=0.;
	286	sigma2Match =0.;
	287	}
	288
	289	fShapesVar[9] = ptMatch;
	290	fShapesVar[10] = ptDMatch;
	291	fShapesVar[11] = massMatch;
	292	fShapesVar[12] = constMatch;
	293	fShapesVar[13] = angulMatch;
	294	fShapesVar[14] = circMatch;
	295	fShapesVar[15] = lesubMatch;
	296	fShapesVar[16] = sigma2Match;
	297
cfa7b4df	298	fTreeObservableTagging->Fill();
b929303d	299	fTreeObservableTagging->SetWeight(kWeight);
cfa7b4df	300	}
	301
	302	}
	303
	304	return kTRUE;
	305	}
	306
	307	//________________________________________________________________________
e5aca8ab	308	Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0) {
cfa7b4df	309	//calc subtracted jet mass
e5aca8ab	310	if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
cfa7b4df	311	return jet->GetSecondOrderSubtracted();
	312	else
	313	return jet->M();
	314	}
	315
	316	//________________________________________________________________________
a68b1989	317	Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
cfa7b4df	318
a68b1989	319	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	320	if (!jet->GetNumberOfTracks())
	321	return 0;
	322	Double_t den=0.;
	323	Double_t num = 0.;
	324	AliVParticle *vp1 = 0x0;
	325	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	326	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	327	Double_t dphi = vp1->Phi()-jet->Phi();
	328	if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
	329	if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
	330	Double_t dr2 = (vp1->Eta()-jet->Eta())(vp1->Eta()-jet->Eta()) + dphidphi;
	331	Double_t dr = TMath::Sqrt(dr2);
	332	num=num+vp1->Pt()*dr;
	333	den=den+vp1->Pt();
	334	}
	335	return num/den;
	336	}
	337
	338	//________________________________________________________________________
a68b1989	339	Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) {
cfa7b4df	340
e5aca8ab	341	if((fJetShapeSub==kDerivSub) && (jetContNb==0))
cfa7b4df	342	return jet->GetSecondOrderSubtractedAngularity();
cfa7b4df	343	else
a68b1989	344	return Angularity(jet, jetContNb);
cfa7b4df	345
	346	}
	347
	348
	349	//________________________________________________________________________
a68b1989	350	Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
cfa7b4df	351
a68b1989	352	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	353	if (!jet->GetNumberOfTracks())
	354	return 0;
	355	Double_t den=0.;
	356	Double_t num = 0.;
	357	AliVParticle *vp1 = 0x0;
	358	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	359	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	360	num=num+vp1->Pt()*vp1->Pt();
	361	den=den+vp1->Pt();
	362	}
	363	return TMath::Sqrt(num)/den;
	364	}
	365
	366	//________________________________________________________________________
a68b1989	367	Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) {
cfa7b4df	368	//calc subtracted jet mass
e5aca8ab	369	if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
cfa7b4df	370	return jet->GetSecondOrderSubtractedpTD();
cfa7b4df	371	else
a68b1989	372	return PTD(jet, jetContNb);
cfa7b4df	373
	374	}
	375
	376	//_____________________________________________________________________________
a68b1989	377	Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
cfa7b4df	378
a68b1989	379	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	380	if (!jet->GetNumberOfTracks())
	381	return 0;
	382	Double_t mxx = 0.;
	383	Double_t myy = 0.;
	384	Double_t mxy = 0.;
	385	int nc = 0;
	386	Double_t sump2 = 0.;
	387	Double_t pxjet=jet->Px();
	388	Double_t pyjet=jet->Py();
	389	Double_t pzjet=jet->Pz();
	390
	391
	392	//2 general normalized vectors perpendicular to the jet
	393	TVector3 ppJ1(pxjet, pyjet, pzjet);
	394	TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
	395	ppJ3.SetMag(1.);
	396	TVector3 ppJ2(-pyjet, pxjet, 0);
	397	ppJ2.SetMag(1.);
	398	AliVParticle *vp1 = 0x0;
	399	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	400	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	401
	402
	403	TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
	404
	405	//local frame
	406	TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
	407	TVector3 pPerp = pp - pLong;
	408	//projection onto the two perpendicular vectors defined above
	409
	410	Float_t ppjX = pPerp.Dot(ppJ2);
	411	Float_t ppjY = pPerp.Dot(ppJ3);
	412	Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
	413	if(ppjT<=0) return 0;
	414
	415	mxx += (ppjX * ppjX / ppjT);
	416	myy += (ppjY * ppjY / ppjT);
	417	mxy += (ppjX * ppjY / ppjT);
	418	nc++;
	419	sump2 += ppjT;}
	420
	421	if(nc<2) return 0;
	422	if(sump2==0) return 0;
	423	// Sphericity Matrix
	424	Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
	425	TMatrixDSym m0(2,ele);
	426
	427	// Find eigenvectors
	428	TMatrixDSymEigen m(m0);
	429	TVectorD eval(2);
	430	TMatrixD evecm = m.GetEigenVectors();
	431	eval = m.GetEigenValues();
	432	// Largest eigenvector
	433	int jev = 0;
	434	// cout<<eval[0]<<" "<<eval[1]<<endl;
	435	if (eval[0] < eval[1]) jev = 1;
	436	TVectorD evec0(2);
	437	// Principle axis
	438	evec0 = TMatrixDColumn(evecm, jev);
	439	Double_t compx=evec0[0];
	440	Double_t compy=evec0[1];
	441	TVector2 evec(compx, compy);
	442	Double_t circ=0;
	443	if(jev==1) circ=2*eval[0];
444	if(jev==0) circ=2*eval[1];
445
446	return circ;
447
448
449
450	}
451
452
453
454
455	//________________________________________________________________________
a68b1989	456	Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) {
cfa7b4df	457	//calc subtracted jet mass
cfa7b4df	458
e5aca8ab	459	if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
cfa7b4df	460	return jet->GetSecondOrderSubtractedCircularity();
cfa7b4df	461	else
a68b1989	462	return Circularity(jet, jetContNb);
cfa7b4df	463
	464	}
	465
	466	//________________________________________________________________________
a68b1989	467	Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
cfa7b4df	468
a68b1989	469	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	470	if (!jet->GetNumberOfTracks())
	471	return 0;
	472	Double_t den=0.;
	473	Double_t num = 0.;
	474	AliVParticle *vp1 = 0x0;
	475	AliVParticle *vp2 = 0x0;
	476	std::vector<int> ordindex;
	477	ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
	478
	479	vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
	480	vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
	481
	482	num=vp1->Pt();
	483	den=vp2->Pt();
	484
	485	return num-den;
	486	}
	487
	488	//________________________________________________________________________
a68b1989	489	Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
cfa7b4df	490	//calc subtracted jet mass
cfa7b4df	491
e5aca8ab	492	if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
cfa7b4df	493	return jet->GetSecondOrderSubtractedLeSub();
cfa7b4df	494	else
a68b1989	495	return LeSub(jet, jetContNb);
cfa7b4df	496
	497	}
	498
	499	//________________________________________________________________________
e5aca8ab	500	Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0) {
cfa7b4df	501	//calc subtracted jet mass
cfa7b4df	502
e5aca8ab	503	if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
cfa7b4df	504	return jet->GetSecondOrderSubtractedConstituent();
	505	else
	506	return jet->GetNumberOfTracks();
	507
	508	}
	509
	510
	511	//______________________________________________________________________________
a68b1989	512	Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
cfa7b4df	513
a68b1989	514	AliJetContainer *jetCont = GetJetContainer(jetContNb);
cfa7b4df	515	if (!jet->GetNumberOfTracks())
	516	return 0;
	517	Double_t mxx = 0.;
	518	Double_t myy = 0.;
	519	Double_t mxy = 0.;
	520	int nc = 0;
	521	Double_t sump2 = 0.;
	522
	523	AliVParticle *vp1 = 0x0;
	524	for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
	525	vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
	526	Double_t ppt=vp1->Pt();
	527	Double_t dphi = vp1->Phi()-jet->Phi();
	528	if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
	529	if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
	530	Double_t deta = vp1->Eta()-jet->Eta();
	531	mxx += pptpptdeta*deta;
	532	myy += pptpptdphi*dphi;
	533	mxy -= pptpptdeta*dphi;
	534	nc++;
	535	sump2 += ppt*ppt;
	536
	537	}
	538	if(nc<2) return 0;
	539	if(sump2==0) return 0;
	540	// Sphericity Matrix
	541	Double_t ele[4] = {mxx , mxy , mxy , myy };
	542	TMatrixDSym m0(2,ele);
	543
	544	// Find eigenvectors
	545	TMatrixDSymEigen m(m0);
	546	TVectorD eval(2);
	547	TMatrixD evecm = m.GetEigenVectors();
	548	eval = m.GetEigenValues();
	549	// Largest eigenvector
	550	int jev = 0;
	551	// cout<<eval[0]<<" "<<eval[1]<<endl;
	552	if (eval[0] < eval[1]) jev = 1;
	553	TVectorD evec0(2);
	554	// Principle axis
	555	evec0 = TMatrixDColumn(evecm, jev);
	556	Double_t compx=evec0[0];
	557	Double_t compy=evec0[1];
	558	TVector2 evec(compx, compy);
	559	Double_t sig=0;
	560	if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
	561	if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
	562
	563	return sig;
	564
	565	}
	566
	567	//________________________________________________________________________
a68b1989	568	Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) {
cfa7b4df	569	//calc subtracted jet mass
cfa7b4df	570
e5aca8ab	571	if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
cfa7b4df	572	return jet->GetSecondOrderSubtractedSigma2();
cfa7b4df	573	else
a68b1989	574	return Sigma2(jet, jetContNb);
cfa7b4df	575
	576	}
	577
	578	//________________________________________________________________________
	579	Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
	580	//
	581	// retrieve event objects
	582	//
	583	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	584	return kFALSE;
	585
	586	return kTRUE;
	587	}
	588
	589	//_______________________________________________________________________
	590	void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *)
	591	{
	592	// Called once at the end of the analysis.
	593
	594	// fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
	595	// if (!fTreeObservableTagging){
	596	// Printf("ERROR: fTreeObservableTagging not available");
	597	// return;
	598	// }
	599
	600	}
	601