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