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