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