PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.cxx

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