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   fPhiJetCorr6(0x0),
  61   fPhiJetCorr7(0x0),
  62   fEtaJetCorr6(0x0),
  63   fEtaJetCorr7(0x0),
  64   fPtJetCorr(0x0),
  65   fPtJet(0x0),
  66   fTreeObservableTagging(0)
  67
  68 {
  69   SetMakeGeneralHistograms(kTRUE);
  70 }
  71
  72 //________________________________________________________________________
  73 AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) :
  74   AliAnalysisTaskEmcalJet(name, kTRUE),
  75   fContainer(0),
  76   fMinFractionShared(0),
  77   fJetShapeType(kData),
  78   fJetShapeSub(kNoSub),
  79   fShapesVar(0),
  80   fPtThreshold(-9999.),
  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)
  89
  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");
 116   Int_t nVar = 17;
 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";
 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";
 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, 2*TMath::Pi(), 50, 0, 2*TMath::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, 2*TMath::Pi(), 50, 0, 2*TMath::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);
 184   Float_t kWeight=1;
 185   if(jetCont) {
 186     jetCont->ResetCurrentID();
 187     while((jet1 = jetCont->GetNextAcceptJet())) {
 188       if (!jet1) continue;
 189       AliEmcalJet* jet2 = 0x0;
 190       fPtJet->Fill(jet1->Pt());
 191
 192       if (!(fJetShapeType == kData)) {
 193         AliPythiaInfo *partonsInfo = 0x0;
 194         if((fJetShapeType == kTrueDet) || (fJetShapeType == kDetEmb)){
 195           AliJetContainer *jetContTrue = GetJetContainer(1);
 196           jet2 = jet1->ClosestJet();
 197           if (!jet2) {
 198             Printf("jet2 not exists, returning");
 199             continue;
 200           }
 201
 202           Double_t fraction = jetCont->GetFractionSharedPt(jet1);
 203           cout<<"hey a jet"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<endl;
 204           if(fraction<fMinFractionShared) continue;
 205           partonsInfo = (AliPythiaInfo*) jetContTrue->GetPythiaInfo();
 206           if(!partonsInfo) return 0;
 207
 208         }
 209         else {
 210           partonsInfo = (AliPythiaInfo*) jetCont->GetPythiaInfo();
 211           jet2=jet1;
 212           if(!partonsInfo) return 0;
 213         }
 214
 215         Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6());
 216         Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
 217         kWeight=partonsInfo->GetPythiaEventWeight();
 218         if (jp1< -1*TMath::Pi()) jp1 = (-2*TMath::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< -1*TMath::Pi()) jp1= (-2*TMath::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
 245       Double_t ptSubtracted = 0;
 246
 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
 254       fShapesVar[1] = ptSubtracted;
 255       fShapesVar[2] = GetJetpTD(jet1,0);
 256       fShapesVar[3] = GetJetMass(jet1,0);
 257       fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
 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);
 269         massMatch=GetJetMass(jet2,kMatched);
 270         constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
 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
 298       fTreeObservableTagging->Fill();
 299       fTreeObservableTagging->SetWeight(kWeight);
 300     }
 301
 302   }
 303
 304   return kTRUE;
 305 }
 306
 307 //________________________________________________________________________
 308 Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0) {
 309   //calc subtracted jet mass
 310   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 311     return jet->GetSecondOrderSubtracted();
 312   else
 313     return jet->M();
 314 }
 315
 316 //________________________________________________________________________
 317 Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
 318
 319   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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()) + dphi*dphi;
 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 //________________________________________________________________________
 339 Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) {
 340
 341   if((fJetShapeSub==kDerivSub) && (jetContNb==0))
 342     return jet->GetSecondOrderSubtractedAngularity();
 343   else
 344     return Angularity(jet, jetContNb);
 345
 346 }
 347
 348
 349 //________________________________________________________________________
 350 Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
 351
 352   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 //________________________________________________________________________
 367 Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) {
 368   //calc subtracted jet mass
 369   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 370     return jet->GetSecondOrderSubtractedpTD();
 371   else
 372     return PTD(jet, jetContNb);
 373
 374 }
 375
 376 //_____________________________________________________________________________
 377 Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
 378
 379   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 //________________________________________________________________________
 456 Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) {
 457   //calc subtracted jet mass
 458
 459   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 460     return jet->GetSecondOrderSubtractedCircularity();
 461   else
 462     return Circularity(jet, jetContNb);
 463
 464 }
 465
 466 //________________________________________________________________________
 467 Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
 468
 469   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 //________________________________________________________________________
 489 Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
 490   //calc subtracted jet mass
 491
 492   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 493     return jet->GetSecondOrderSubtractedLeSub();
 494   else
 495     return LeSub(jet, jetContNb);
 496
 497 }
 498
 499 //________________________________________________________________________
 500 Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0) {
 501   //calc subtracted jet mass
 502
 503   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 504     return jet->GetSecondOrderSubtractedConstituent();
 505   else
 506     return jet->GetNumberOfTracks();
 507
 508 }
 509
 510
 511 //______________________________________________________________________________
 512 Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
 513
 514   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 += ppt*ppt*deta*deta;
 532        myy += ppt*ppt*dphi*dphi;
 533        mxy -= ppt*ppt*deta*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 //________________________________________________________________________
 568 Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) {
 569   //calc subtracted jet mass
 570
 571   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 572     return jet->GetSecondOrderSubtractedSigma2();
 573   else
 574     return Sigma2(jet, jetContNb);
 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