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