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 "AliStackPartonInfo.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
 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         AliStackPartonInfo *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 = (AliStackPartonInfo*) jetContTrue->GetStackPartonInfo();
 206           if(!partonsInfo) return 0;
 207         }
 208         else {
 209           partonsInfo = (AliStackPartonInfo*) jetCont->GetStackPartonInfo();
 210           jet2=jet1;
 211           if(!partonsInfo) return 0;
 212         }
 213
 214         Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6());
 215         Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
 216
 217         if (jp1< -1*TMath::Pi()) jp1 = (-2*TMath::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< -1*TMath::Pi()) jp1= (-2*TMath::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
 244       Double_t ptSubtracted = 0;
 245
 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
 253       fShapesVar[1] = ptSubtracted;
 254       fShapesVar[2] = GetJetpTD(jet1,0);
 255       fShapesVar[3] = GetJetMass(jet1,0);
 256       fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
 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,kMatched);
 269         constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
 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
 297       fTreeObservableTagging->Fill();
 298
 299     }
 300
 301   }
 302
 303   return kTRUE;
 304 }
 305
 306 //________________________________________________________________________
 307 Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0) {
 308   //calc subtracted jet mass
 309   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 310     return jet->GetSecondOrderSubtracted();
 311   else
 312     return jet->M();
 313 }
 314
 315 //________________________________________________________________________
 316 Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
 317
 318   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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()) + dphi*dphi;
 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 //________________________________________________________________________
 338 Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) {
 339
 340   if((fJetShapeSub==kDerivSub) && (jetContNb==0))
 341     return jet->GetSecondOrderSubtractedAngularity();
 342   else
 343     return Angularity(jet, jetContNb);
 344
 345 }
 346
 347
 348 //________________________________________________________________________
 349 Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
 350
 351   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 //________________________________________________________________________
 366 Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) {
 367   //calc subtracted jet mass
 368   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 369     return jet->GetSecondOrderSubtractedpTD();
 370   else
 371     return PTD(jet, jetContNb);
 372
 373 }
 374
 375 //_____________________________________________________________________________
 376 Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
 377
 378   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 //________________________________________________________________________
 455 Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) {
 456   //calc subtracted jet mass
 457
 458   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 459     return jet->GetSecondOrderSubtractedCircularity();
 460   else
 461     return Circularity(jet, jetContNb);
 462
 463 }
 464
 465 //________________________________________________________________________
 466 Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
 467
 468   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 //________________________________________________________________________
 488 Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
 489   //calc subtracted jet mass
 490
 491   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 492     return jet->GetSecondOrderSubtractedLeSub();
 493   else
 494     return LeSub(jet, jetContNb);
 495
 496 }
 497
 498 //________________________________________________________________________
 499 Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0) {
 500   //calc subtracted jet mass
 501
 502   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 503     return jet->GetSecondOrderSubtractedConstituent();
 504   else
 505     return jet->GetNumberOfTracks();
 506
 507 }
 508
 509
 510 //______________________________________________________________________________
 511 Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
 512
 513   AliJetContainer *jetCont = GetJetContainer(jetContNb);
 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 += ppt*ppt*deta*deta;
 531        myy += ppt*ppt*dphi*dphi;
 532        mxy -= ppt*ppt*deta*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 //________________________________________________________________________
 567 Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) {
 568   //calc subtracted jet mass
 569
 570   if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
 571     return jet->GetSecondOrderSubtractedSigma2();
 572   else
 573     return Sigma2(jet, jetContNb);
 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