AliPhysics  a9863a5 (a9863a5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalQGTagging.cxx
Go to the documentation of this file.
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 "AliEmcalPythiaInfo.h"
38 #include "TRandom3.h"
39 
40 
41 
42 #include "AliAODEvent.h"
44 
45 using std::cout;
46 using std::endl;
47 
49 
50 //________________________________________________________________________
53  fContainer(0),
54  fMinFractionShared(0),
55  fJetShapeType(kData),
56  fJetShapeSub(kNoSub),
57  fJetSelection(kInclusive),
58  fShapesVar(0),
59  fPtThreshold(-9999.),
60  fRMatching(0.2),
61  fSelectedShapes(0),
62  fminpTTrig(20.),
63  fmaxpTTrig(50.),
64  fangWindowRecoil(0.6),
65  fSemigoodCorrect(0),
66  fHolePos(0),
67  fHoleWidth(0),
68  fCentSelectOn(kTRUE),
69  fCentMin(0),
70  fCentMax(10),
71  fOneConstSelectOn(kFALSE),
72  fDerivSubtrOrder(0),
73  fh2ResponseUW(0x0),
74  fh2ResponseW(0x0),
75  fPhiJetCorr6(0x0),
76  fPhiJetCorr7(0x0),
77  fEtaJetCorr6(0x0),
78  fEtaJetCorr7(0x0),
79  fPtJetCorr(0x0),
80  fPtJet(0x0),
81  fhpTjetpT(0x0),
82  fhPt(0x0),
83  fhPhi(0x0),
84  fNbOfConstvspT(0x0),
85  fTreeObservableTagging(0)
86 
87 {
88  SetMakeGeneralHistograms(kTRUE);
89 }
90 
91 //________________________________________________________________________
93  AliAnalysisTaskEmcalJet(name, kTRUE),
94  fContainer(0),
95  fMinFractionShared(0),
96  fJetShapeType(kData),
97  fJetShapeSub(kNoSub),
98  fJetSelection(kInclusive),
99  fShapesVar(0),
100  fPtThreshold(-9999.),
101  fRMatching(0.2),
102  fSelectedShapes(0),
103  fminpTTrig(20.),
104  fmaxpTTrig(50.),
105  fangWindowRecoil(0.6),
106  fSemigoodCorrect(0),
107  fHolePos(0),
108  fHoleWidth(0),
109  fCentSelectOn(kTRUE),
110  fCentMin(0),
111  fCentMax(10),
112  fOneConstSelectOn(kFALSE),
113  fDerivSubtrOrder(0),
114  fh2ResponseUW(0x0),
115  fh2ResponseW(0x0),
116  fPhiJetCorr6(0x0),
117  fPhiJetCorr7(0x0),
118  fEtaJetCorr6(0x0),
119  fEtaJetCorr7(0x0),
120  fPtJetCorr(0x0),
121  fPtJet(0x0),
122  fhpTjetpT(0x0),
123  fhPt(0x0),
124  fhPhi(0x0),
125  fNbOfConstvspT(0x0),
126  fTreeObservableTagging(0)
127 
128 {
129  // Standard constructor.
130 
132 
133  DefineOutput(1, TTree::Class());
134 
135 }
136 
137 //________________________________________________________________________
139 {
140  // Destructor.
141 }
142 
143 //________________________________________________________________________
145 {
146  // Create user output.
147 
149 
150  Bool_t oldStatus = TH1::AddDirectoryStatus();
151  TH1::AddDirectory(kFALSE);
152 
153  fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
154  const Int_t nVar = 14;
155  const Int_t nVarless = 10;
156 
157  if(fSelectedShapes==0){
158 
159  fShapesVar = new Float_t [nVar];
160  TString *fShapesVarNames = new TString [nVar];
161 
162  fShapesVarNames[0] = "partonCode";
163  fShapesVarNames[1] = "ptJet";
164  fShapesVarNames[2] = "ptDJet";
165  fShapesVarNames[3] = "mJet";
166  // fShapesVarNames[4] = "nbOfConst";
167  fShapesVarNames[4] = "angularity";
168  fShapesVarNames[5] = "circularity";
169  fShapesVarNames[6] = "lesub";
170  //fShapesVarNames[6] = "sigma2";
171 
172  fShapesVarNames[7] = "ptJetMatch";
173  fShapesVarNames[8] = "ptDJetMatch";
174  fShapesVarNames[9] = "mJetMatch";
175  // fShapesVarNames[12] = "nbOfConstMatch";
176  fShapesVarNames[10] = "angularityMatch";
177  fShapesVarNames[11] = "circularityMatch";
178  fShapesVarNames[12] = "lesubMatch";
179  //fShapesVarNames[12] = "sigma2Match";
180  fShapesVarNames[13]="weightPythia";
181 
182  for(Int_t ivar=0; ivar < nVar; ivar++){
183  cout<<"looping over variables"<<endl;
184  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
185 
186 
187 }
188 
189 
190  if(fSelectedShapes==1){
191  fShapesVar = new Float_t [nVarless];
192 
193  TString *fShapesVarNames = new TString [nVarless];
194  fShapesVarNames[0] = "partonCode";
195  fShapesVarNames[1] = "ptJet";
196  fShapesVarNames[2] = "ptDJet";
197  fShapesVarNames[3] = "angularity";
198  fShapesVarNames[4] = "lesub";
199  fShapesVarNames[5] = "ptJetMatch";
200  fShapesVarNames[6] = "ptDJetMatch";
201  fShapesVarNames[7] = "angularityMatch";
202  fShapesVarNames[8] = "lesubMatch";
203  fShapesVarNames[9]="weightPythia";
204 
205  for(Int_t ivar=0; ivar < nVarless; ivar++){
206  cout<<"looping over variables"<<endl;
207  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
208 
209 }
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
222 
223 
224 
225  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
226  fOutput->Add(fh2ResponseUW);
227  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
228  fOutput->Add(fh2ResponseW);
229  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
230  fOutput->Add(fPhiJetCorr6);
231  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
232  fOutput->Add(fEtaJetCorr6);
233 
234  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
235  fOutput->Add(fPhiJetCorr7);
236  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
237  fOutput->Add(fEtaJetCorr7);
238 
239  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
240  fOutput->Add(fPtJetCorr);
241  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
242  fOutput->Add(fPtJet);
243 
244  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
245  fOutput->Add(fhpTjetpT);
246  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
247  fOutput->Add(fhPt);
248  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
249  fOutput->Add(fhPhi);
250 
251  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
252  fOutput->Add(fNbOfConstvspT);
253 
255  TH1::AddDirectory(oldStatus);
256  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
257 
258 }
259 
260 //________________________________________________________________________
262 {
263  // Run analysis code here, if needed. It will be executed before FillHistograms().
264 
265  return kTRUE;
266 }
267 
268 //________________________________________________________________________
270 {
271  // Fill histograms.
272  //cout<<"base container"<<endl;
273  AliEmcalJet* jet1 = NULL;
274  AliJetContainer *jetCont = GetJetContainer(0);
275 
276  Float_t kWeight=1;
277  if (fCentSelectOn)
278  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
279 
280  AliAODTrack *triggerHadron = 0x0;
281 
282  if (fJetSelection == kRecoil) {
283  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
284  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
285 
286 
287  if (triggerHadronLabel==-99999) {
288  //Printf ("Trigger Hadron not found, return");
289  return 0;}
290 
291 
293  TClonesArray *trackArrayAn = partContAn->GetArray();
294  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
295 
296  if (!triggerHadron) {
297  //Printf("No Trigger hadron with the found label!!");
298  return 0;
299  }
300 
301  if(fSemigoodCorrect){
302  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
303  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
304  return 0;}
305  }
306 
307  fhPt->Fill(triggerHadron->Pt());
308 
309  }
310 
311  if(jetCont) {
312  jetCont->ResetCurrentID();
313  while((jet1 = jetCont->GetNextAcceptJet())) {
314  if (!jet1) continue;
315  AliEmcalJet* jet2 = 0x0;
316  AliEmcalJet* jet3 = 0x0;
317  fPtJet->Fill(jet1->Pt());
318  AliEmcalJet *jetUS = NULL;
319  Int_t ifound=0, jfound=0;
320  Int_t ilab=-1, jlab=-1;
321 
323  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
324  if(TMath::Abs(disthole)<fHoleWidth){
325  continue;}
326  }
327 
328  Float_t dphiRecoil = 0.;
329  if (fJetSelection == kRecoil){
330  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
331  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
332  // Printf("Recoil jets back to back not found! continuing");
333  continue;
334  }
335 
336  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
337  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
338  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
339 
340  }
341 
342 
343  fShapesVar[0] = 0.;
345  AliJetContainer *jetContTrue = GetJetContainer(1);
346  AliJetContainer *jetContUS = GetJetContainer(2);
347 
348  if(fJetShapeSub==kConstSub){
349 
350  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
351  jetUS = jetContUS->GetJet(i);
352  if(jetUS->GetLabel()==jet1->GetLabel()) {
353  ifound++;
354  if(ifound==1) ilab = i;
355  }
356  }
357  if(ilab==-1) continue;
358  jetUS=jetContUS->GetJet(ilab);
359  jet2=jetUS->ClosestJet();
360  }
361 
362  if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet();
363  if (!jet2) {
364  Printf("jet2 does not exist, returning");
365  continue;
366  }
367 
368  AliJetContainer *jetContPart=GetJetContainer(3);
369  jet3=jet2->ClosestJet();
370 
371  if(!jet3){
372  Printf("jet3 does not exist, returning");
373  continue;
374  }
375  cout<<"jet 3 exists"<<jet3->Pt()<<endl;
376 
377 
378  fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt());
379 
380  Double_t fraction=0;
381  if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1);
382  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
383  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
384 
385  if(fraction<fMinFractionShared) continue;
386  //InputEvent()->Print();
387 
388  }
389 
390 
391  if (fJetShapeType == kPythiaDef){
392 
393  AliJetContainer *jetContTrue = GetJetContainer(1);
394  AliJetContainer *jetContUS = GetJetContainer(2);
395  AliJetContainer *jetContPart = GetJetContainer(3);
396 
397  if(fJetShapeSub==kConstSub){
398 
399  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
400  jetUS = jetContUS->GetJet(i);
401  if(jetUS->GetLabel()==jet1->GetLabel()) {
402  ifound++;
403  if(ifound==1) ilab = i;
404  }
405  }
406  if(ilab==-1) continue;
407  jetUS=jetContUS->GetJet(ilab);
408  jet2=jetUS->ClosestJet();
409 
410  if (!jet2) {
411  Printf("jet2 does not exist, returning");
412  continue;
413  }
414 
415  for(Int_t j=0; j<jetContPart->GetNJets(); j++) {
416 
417  jet3 = jetContPart->GetJet(j);
418  if(!jet3) continue;
419  if(jet3->GetLabel()==jet2->GetLabel()) {
420  jfound++;
421  if(jfound==1) jlab = j;
422  }
423  }
424  if(jlab==-1) continue;
425  jet3=jetContPart->GetJet(jlab);
426  if(!jet3){
427  Printf("jet3 does not exist, returning");
428  continue;
429  }
430  }
431  if(!(fJetShapeSub==kConstSub)) jet3 = jet1->ClosestJet();
432  if (!jet3) {
433  Printf("jet3 does not exist, returning");
434  continue;
435  }
436 
437 
438  fh2ResponseUW->Fill(jet1->Pt(),jet3->Pt());
439 
440 
441  }
442 
443 
444  if (fJetShapeType == kGenOnTheFly){
445  const AliEmcalPythiaInfo *partonsInfo = 0x0;
446  partonsInfo = GetPythiaInfo();
447  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
448  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
449  kWeight=partonsInfo->GetPythiaEventWeight();
450  fh2ResponseW->Fill(jet1->Pt(),jet1->Pt(),kWeight);
451 
452  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
453  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
454  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
455  if(dRp1 < fRMatching) {
456  fShapesVar[0] = partonsInfo->GetPartonFlag6();
457  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
458  }
459  else {
460  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
461  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
462  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
463  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
464  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
465  if(dRp1 < fRMatching) {
466  fShapesVar[0] = partonsInfo->GetPartonFlag7();
467  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
468  }
469  else fShapesVar[0]=0;
470  }
471  }
472 
473 
474 
475 
476  Double_t ptSubtracted = 0;
477  if (fJetShapeSub==kConstSub) ptSubtracted= jet1->Pt();
478 
479  else if (fJetShapeSub==kDerivSub) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
480 
481  else if (fJetShapeSub==kNoSub) {
482  if ((fJetShapeType==kData) || (fJetShapeType==kDetEmbPartPythia)) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
483  else if ((fJetShapeType==kPythiaDef) || (fJetShapeType==kMCTrue) || (fJetShapeType==kGenOnTheFly)) ptSubtracted= jet1->Pt();
484  }
485  if (ptSubtracted < fPtThreshold) continue;
486 
487  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
488  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
489 
490 
491 
492  if(fSelectedShapes==0){
493  fShapesVar[1] = ptSubtracted;
494  fShapesVar[2] = GetJetpTD(jet1,0);
495  fShapesVar[3] = GetJetMass(jet1,0);
496  fShapesVar[4] = GetJetAngularity(jet1,0);
497  fShapesVar[5] = GetJetCircularity(jet1,0);
498  fShapesVar[6] = GetJetLeSub(jet1,0);}
499 
500  if(fSelectedShapes==1){
501  fShapesVar[1] = ptSubtracted;
502  fShapesVar[2] = GetJetpTD(jet1,0);
503  fShapesVar[3] = GetJetAngularity(jet1,0);
504  fShapesVar[4] = GetJetLeSub(jet1,0);}
505 
506 
507  Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0.;
508  Int_t kMatched = 0;
509 
510  if (fJetShapeType==kPythiaDef) {
511  kMatched =1;
512  if(fJetShapeSub==kConstSub) kMatched = 3;
513 
514  ptMatch=jet3->Pt();
515  ptDMatch=GetJetpTD(jet3, kMatched);
516  if(fSelectedShapes==0) massMatch=GetJetMass(jet3,kMatched);
517  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
518  angulMatch=GetJetAngularity(jet3, kMatched);
519  if(fSelectedShapes==0) circMatch=GetJetCircularity(jet3, kMatched);
520  lesubMatch=GetJetLeSub(jet3, kMatched);
521  //sigma2Match = GetSigma2(jet2, kMatched);
522  }
523 
525  if(fJetShapeSub==kConstSub) kMatched = 3;
526  if(fJetShapeSub==kDerivSub) kMatched = 2;
527  ptMatch=jet3->Pt();
528  ptDMatch=GetJetpTD(jet3, kMatched);
529  if(fSelectedShapes==0) massMatch=GetJetMass(jet3,kMatched);
530  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
531  angulMatch=GetJetAngularity(jet3, kMatched);
532  if(fSelectedShapes==0) circMatch=GetJetCircularity(jet3, kMatched);
533  lesubMatch=GetJetLeSub(jet3, kMatched);
534  //sigma2Match = GetSigma2(jet3, kMatched)
535 
536  }
537 
538 
539 
541  kMatched = 0;
542  ptMatch=0.;
543  ptDMatch=0.;
544  massMatch=0.;
545  //constMatch=0.;
546  angulMatch=0.;
547  circMatch=0.;
548  lesubMatch=0.;
549  //sigma2Match =0.;
550 
551  }
552 
553  if(fSelectedShapes==0){
554 
555  fShapesVar[7] = ptMatch;
556  fShapesVar[8] = ptDMatch;
557  fShapesVar[9] = massMatch;
558  fShapesVar[10] = angulMatch;
559  fShapesVar[11] = circMatch;
560  fShapesVar[12] = lesubMatch;
561  fShapesVar[13] = kWeight;
562 
563  }
564 
565 
566  if(fSelectedShapes==1){
567 
568  fShapesVar[5] = ptMatch;
569  fShapesVar[6] = ptDMatch;
570  fShapesVar[7] = angulMatch;
571  fShapesVar[8] = lesubMatch;
572  fShapesVar[9] = kWeight;}
573 
574 
575  fTreeObservableTagging->Fill();
576 
577 
578 
579 
580 
581 
582  }
583 
584  }
585 
586  return kTRUE;
587 }
588 
589 //________________________________________________________________________
591  //calc subtracted jet mass
592  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
594  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
595  else
596  return jet->M();
597 }
598 
599 //________________________________________________________________________
600 Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
601 
602  AliJetContainer *jetCont = GetJetContainer(jetContNb);
603  if (!jet->GetNumberOfTracks())
604  return 0;
605  Double_t den=0.;
606  Double_t num = 0.;
607  AliVParticle *vp1 = 0x0;
608  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
609  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
610 
611  if (!vp1){
612  Printf("AliVParticle associated to constituent not found");
613  continue;
614  }
615 
616  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
617  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
618  Double_t dr = TMath::Sqrt(dr2);
619  num=num+vp1->Pt()*dr;
620  den=den+vp1->Pt();
621  }
622  return num/den;
623 }
624 
625 //________________________________________________________________________
627 
628  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
631  else
632  return Angularity(jet, jetContNb);
633 
634 }
635 
636 
637 //________________________________________________________________________
638 Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
639 
640  AliJetContainer *jetCont = GetJetContainer(jetContNb);
641  if (!jet->GetNumberOfTracks())
642  return 0;
643  Double_t den=0.;
644  Double_t num = 0.;
645  AliVParticle *vp1 = 0x0;
646  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
647  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
648 
649  if (!vp1){
650  Printf("AliVParticle associated to constituent not found");
651  continue;
652  }
653 
654  num=num+vp1->Pt()*vp1->Pt();
655  den=den+vp1->Pt();
656  }
657  return TMath::Sqrt(num)/den;
658 }
659 
660 //________________________________________________________________________
661 Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0){
662  //calc subtracted jet mass
663  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
665  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
666  else
667  return PTD(jet, jetContNb);
668 
669 }
670 
671 //_____________________________________________________________________________
672 Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
673 
674  AliJetContainer *jetCont = GetJetContainer(jetContNb);
675  if (!jet->GetNumberOfTracks())
676  return 0;
677  Double_t mxx = 0.;
678  Double_t myy = 0.;
679  Double_t mxy = 0.;
680  int nc = 0;
681  Double_t sump2 = 0.;
682  Double_t pxjet=jet->Px();
683  Double_t pyjet=jet->Py();
684  Double_t pzjet=jet->Pz();
685 
686 
687  //2 general normalized vectors perpendicular to the jet
688  TVector3 ppJ1(pxjet, pyjet, pzjet);
689  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
690  ppJ3.SetMag(1.);
691  TVector3 ppJ2(-pyjet, pxjet, 0);
692  ppJ2.SetMag(1.);
693  AliVParticle *vp1 = 0x0;
694  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
695  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
696 
697  if (!vp1){
698  Printf("AliVParticle associated to constituent not found");
699  continue;
700  }
701 
702  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
703 
704  //local frame
705  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
706  TVector3 pPerp = pp - pLong;
707  //projection onto the two perpendicular vectors defined above
708 
709  Float_t ppjX = pPerp.Dot(ppJ2);
710  Float_t ppjY = pPerp.Dot(ppJ3);
711  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
712  if(ppjT<=0) return 0;
713 
714  mxx += (ppjX * ppjX / ppjT);
715  myy += (ppjY * ppjY / ppjT);
716  mxy += (ppjX * ppjY / ppjT);
717  nc++;
718  sump2 += ppjT;}
719 
720  if(nc<2) return 0;
721  if(sump2==0) return 0;
722  // Sphericity Matrix
723  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
724  TMatrixDSym m0(2,ele);
725 
726  // Find eigenvectors
727  TMatrixDSymEigen m(m0);
728  TVectorD eval(2);
729  TMatrixD evecm = m.GetEigenVectors();
730  eval = m.GetEigenValues();
731  // Largest eigenvector
732  int jev = 0;
733  // cout<<eval[0]<<" "<<eval[1]<<endl;
734  if (eval[0] < eval[1]) jev = 1;
735  TVectorD evec0(2);
736  // Principle axis
737  evec0 = TMatrixDColumn(evecm, jev);
738  Double_t compx=evec0[0];
739  Double_t compy=evec0[1];
740  TVector2 evec(compx, compy);
741  Double_t circ=0;
742  if(jev==1) circ=2*eval[0];
743  if(jev==0) circ=2*eval[1];
744 
745  return circ;
746 
747 
748 
749 }
750 
751 
752 
753 
754 //________________________________________________________________________
756  //calc subtracted jet mass
757 
758  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
761  else
762  return Circularity(jet, jetContNb);
763 
764 }
765 
766 //________________________________________________________________________
767 Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
768 
769  AliJetContainer *jetCont = GetJetContainer(jetContNb);
770  if (!jet->GetNumberOfTracks())
771  return 0;
772  Double_t den=0.;
773  Double_t num = 0.;
774  AliVParticle *vp1 = 0x0;
775  AliVParticle *vp2 = 0x0;
776  std::vector<int> ordindex;
777  ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
778  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
779  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
780 
781  if(ordindex.size()<2) return -1;
782 
783  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
784  if (!vp1){
785  Printf("AliVParticle associated to Leading constituent not found");
786  return -1;
787  }
788 
789  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
790  if (!vp2){
791  Printf("AliVParticle associated to Subleading constituent not found");
792  return -1;
793  }
794 
795 
796  num=vp1->Pt();
797  den=vp2->Pt();
798  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
799 
800 return num-den;
801 }
802 
803 //________________________________________________________________________
804 Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
805  //calc subtracted jet mass
806 
807  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
810  else
811  return LeSub(jet, jetContNb);
812 
813 }
814 
815 //________________________________________________________________________
817  //calc subtracted jet mass
818 
819  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
822  else
823  return jet->GetNumberOfTracks();
824 
825  }
826 
827 
828 //______________________________________________________________________________
829 Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
830 
831  AliJetContainer *jetCont = GetJetContainer(jetContNb);
832  if (!jet->GetNumberOfTracks())
833  return 0;
834  Double_t mxx = 0.;
835  Double_t myy = 0.;
836  Double_t mxy = 0.;
837  int nc = 0;
838  Double_t sump2 = 0.;
839 
840  AliVParticle *vp1 = 0x0;
841  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
842  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
843 
844  if (!vp1){
845  Printf("AliVParticle associated to constituent not found");
846  continue;
847  }
848 
849  Double_t ppt=vp1->Pt();
850  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
851 
852  Double_t deta = vp1->Eta()-jet->Eta();
853  mxx += ppt*ppt*deta*deta;
854  myy += ppt*ppt*dphi*dphi;
855  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
856  nc++;
857  sump2 += ppt*ppt;
858 
859  }
860  if(nc<2) return 0;
861  if(sump2==0) return 0;
862  // Sphericity Matrix
863  Double_t ele[4] = {mxx , mxy , mxy , myy };
864  TMatrixDSym m0(2,ele);
865 
866  // Find eigenvectors
867  TMatrixDSymEigen m(m0);
868  TVectorD eval(2);
869  TMatrixD evecm = m.GetEigenVectors();
870  eval = m.GetEigenValues();
871  // Largest eigenvector
872  int jev = 0;
873  // cout<<eval[0]<<" "<<eval[1]<<endl;
874  if (eval[0] < eval[1]) jev = 1;
875  TVectorD evec0(2);
876  // Principle axis
877  evec0 = TMatrixDColumn(evecm, jev);
878  Double_t compx=evec0[0];
879  Double_t compy=evec0[1];
880  TVector2 evec(compx, compy);
881  Double_t sig=0;
882  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
883  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
884 
885  return sig;
886 
887 }
888 
889 //________________________________________________________________________
890 Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0){
891  //calc subtracted jet mass
892 
893  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
896  else
897  return Sigma2(jet, jetContNb);
898 
899 }
900 
901 //________________________________________________________________________
902 Int_t AliAnalysisTaskEmcalQGTagging::SelectTrigger(Float_t minpT, Float_t maxpT){
903 
905  TClonesArray *tracksArray = partCont->GetArray();
906 
907  if(!partCont || !tracksArray) return -99999;
908  AliAODTrack *track = 0x0;
909  AliEmcalParticle *emcPart = 0x0;
910 
911 
912  TList *trackList = new TList();
913  Int_t triggers[100];
914  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
915  Int_t iTT = 0;
916 
917  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
918 
919 
920  if (fJetShapeSub == kConstSub){
921  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
922  if (!emcPart) continue;
923  if(TMath::Abs(emcPart->Eta())>0.9) continue;
924  if (emcPart->Pt()<0.15) continue;
925 
926  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
927  trackList->Add(emcPart);
928  triggers[iTT] = iTrack;
929  iTT++;
930  }
931  }
932  else{
933  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
934  if (!track) continue;
935  if(TMath::Abs(track->Eta())>0.9) continue;
936  if (track->Pt()<0.15) continue;
937  if (!(track->TestFilterBit(768))) continue;
938 
939  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
940  trackList->Add(track);
941  triggers[iTT] = iTrack;
942  iTT++;
943 
944  }
945  }
946  }
947 
948  if (iTT == 0) return -99999;
949  Int_t nbRn = 0, index = 0 ;
950  TRandom3* random = new TRandom3(0);
951  nbRn = random->Integer(iTT);
952 
953  index = triggers[nbRn];
954  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
955  return index;
956 
957 }
958 
959 //__________________________________________________________________________________
960 Double_t AliAnalysisTaskEmcalQGTagging::RelativePhi(Double_t mphi,Double_t vphi){
961 
962  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
963  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
964  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
965  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
966  double dphi = mphi-vphi;
967  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
968  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
969  return dphi;//dphi in [-Pi, Pi]
970 }
971 
972 
973 //________________________________________________________________________
975  //
976  // retrieve event objects
977  //
979  return kFALSE;
980 
981  return kTRUE;
982 }
983 
984 //_______________________________________________________________________
986 {
987  // Called once at the end of the analysis.
988 
989  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
990  // if (!fTreeObservableTagging){
991  // Printf("ERROR: fTreeObservableTagging not available");
992  // return;
993  // }
994 
995 }
996 
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Double_t GetFirstOrderSubtractedAngularity() const
Double_t Area() const
Definition: AliEmcalJet.h:97
Double_t GetSecondOrderSubtractedSigma2() const
Float_t GetPartonEta6() const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:193
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:88
Double_t Py() const
Definition: AliEmcalJet.h:74
Double_t Phi() const
Definition: AliEmcalJet.h:84
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPythiaEventWeight() const
Int_t GetLabel() const
Definition: AliEmcalJet.h:91
Float_t GetPartonPhi7() const
TList * fOutput
!output list
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
Int_t GetPartonFlag7() const
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:106
TString kData
Declare data MC or deltaAOD.
Double_t Px() const
Definition: AliEmcalJet.h:73
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb)
Double_t Eta() const
Double_t Pt() const
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
std::vector< int > SortConstituentsPt(TClonesArray *tracks) const
Int_t GetNJets() const
Double_t GetSecondOrderSubtractedAngularity() const
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb)
Double_t fCent
!event centrality
TClonesArray * GetArray() const
AliEmcalJet * GetNextAcceptJet()
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb)
Double_t Pt() const
Definition: AliEmcalJet.h:76
Double_t GetRhoVal(Int_t i=0) const
Double_t GetFirstOrderSubtractedCircularity() const
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:127
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb)
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Double_t GetSecondOrderSubtractedCircularity() const
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const AliEmcalPythiaInfo * GetPythiaInfo() const
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
ClassImp(AliAnalysisTaskEmcalQGTagging) AliAnalysisTaskEmcalQGTagging
Double_t Pz() const
Definition: AliEmcalJet.h:75
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Declaration of class AliEmcalPythiaInfo.
Double_t GetFirstOrderSubtractedSigma2() const
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb)
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:227
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPartonPt7() const
Double_t M() const
Definition: AliEmcalJet.h:87
void ResetCurrentID(Int_t i=-1)
Container for jet within the EMCAL jet framework.
Float_t GetPartonPt6() const
AliEmcalJet * GetJet(Int_t i) const