AliPhysics  97344c9 (97344c9)
 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<17;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<17;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 
202  TH1::AddDirectory(oldStatus);
203  const Int_t nVar = 18;
204  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
205  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
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  fShapesVarNames[14]="ntrksEvt";
230  fShapesVarNames[15]="rhoVal";
231  fShapesVarNames[16]="rhoMassVal";
232  fShapesVarNames[17]="ptUnsub";
233 
234  for(Int_t ivar=0; ivar < nVar; ivar++){
235  cout<<"looping over variables"<<endl;
236  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
237 
238  PostData(1,fOutput);
239  PostData(2,fTreeObservableTagging);
240 
241  delete [] fShapesVarNames;
242 }
243 
244 //________________________________________________________________________
246 {
247  // Run analysis code here, if needed. It will be executed before FillHistograms().
248 
249  return kTRUE;
250 }
251 
252 //________________________________________________________________________
254 {
255  // Fill histograms.
256  //cout<<"base container"<<endl;
257  AliEmcalJet* jet1 = NULL;
258  AliJetContainer *jetCont = GetJetContainer(0);
259 
260  Float_t kWeight=1;
261  if (fCentSelectOn)
262  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
263 
264  AliAODTrack *triggerHadron = 0x0;
265 
266  if (fJetSelection == kRecoil) {
267  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
268  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
269 
270 
271  if (triggerHadronLabel==-99999) {
272  //Printf ("Trigger Hadron not found, return");
273  return 0;}
274 
275 
277  TClonesArray *trackArrayAn = partContAn->GetArray();
278  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
279 
280  if (!triggerHadron) {
281  //Printf("No Trigger hadron with the found label!!");
282  return 0;
283  }
284 
285  if(fSemigoodCorrect){
286  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
287  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
288  return 0;}
289  }
290 
291  fhPt->Fill(triggerHadron->Pt());
292 
293  }
294 
295 
297  TClonesArray *trackArrayAn = partContAn->GetArray();
298  Int_t ntracksEvt = trackArrayAn->GetEntriesFast();
299 
300  Float_t rhoVal=0, rhoMassVal = 0.;
301  if(jetCont) {
302 
303  jetCont->ResetCurrentID();
305  //rho
306  AliRhoParameter* rhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoSparseR020"));
307  if (!rhoParam) {
308  Printf("%s: Could not retrieve rho %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoName().Data());
309  } else rhoVal = rhoParam->GetVal();
310  //rhom
311  AliRhoParameter* rhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoMassSparseR020"));
312  if (!rhomParam) {
313  Printf("%s: Could not retrieve rho_m %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoMassName().Data());
314  } else rhoMassVal = rhomParam->GetVal();
315  }
316 
317  while((jet1 = jetCont->GetNextAcceptJet())) {
318  if (!jet1) continue;
319  AliEmcalJet* jet2 = 0x0;
320  AliEmcalJet* jet3 = 0x0;
321  fPtJet->Fill(jet1->Pt());
322  AliEmcalJet *jetUS = NULL;
323  Int_t ifound=0, jfound=0;
324  Int_t ilab=-1, jlab=-1;
325 
327  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
328  if(TMath::Abs(disthole)<fHoleWidth){
329  continue;}
330  }
331 
332  Float_t dphiRecoil = 0.;
333  if (fJetSelection == kRecoil){
334  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
335  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
336  // Printf("Recoil jets back to back not found! continuing");
337  continue;
338  }
339 
340  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
341  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
342  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
343 
344  }
345 
346 
347  fShapesVar[0] = 0.;
349  AliJetContainer *jetContTrue = GetJetContainer(1);
350  AliJetContainer *jetContUS = GetJetContainer(2);
351 
352  if(fJetShapeSub==kConstSub){
353  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
354  jetUS = jetContUS->GetJet(i);
355  if(jetUS->GetLabel()==jet1->GetLabel()) {
356  ifound++;
357  if(ifound==1) ilab = i;
358  }
359  }
360  if(ilab==-1) continue;
361  jetUS=jetContUS->GetJet(ilab);
362  jet2=jetUS->ClosestJet();
363  }
364 
365  if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet();
366  if (!jet2) {
367  Printf("jet2 does not exist, returning");
368  continue;
369  }
370 
371  AliJetContainer *jetContPart=GetJetContainer(3);
372  jet3=jet2->ClosestJet();
373 
374  if(!jet3){
375  Printf("jet3 does not exist, returning");
376  continue;
377  }
378  cout<<"jet 3 exists"<<jet3->Pt()<<endl;
379 
380 
381  fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt());
382 
383  Double_t fraction=0;
384  if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1);
385  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
386  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
387 
388  if(fraction<fMinFractionShared) continue;
389  //InputEvent()->Print();
390 
391  }
392 
393 
394 
395  if (fJetShapeType == kPythiaDef){
396 
397  AliJetContainer *jetContTrue = GetJetContainer(1);
398  AliJetContainer *jetContUS = GetJetContainer(2);
399  AliJetContainer *jetContPart = GetJetContainer(3);
400 
401  if(fJetShapeSub==kConstSub){
402 
403  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
404  jetUS = jetContUS->GetJet(i);
405  if(jetUS->GetLabel()==jet1->GetLabel()) {
406  ifound++;
407  if(ifound==1) ilab = i;
408  }
409  }
410  if(ilab==-1) continue;
411  jetUS=jetContUS->GetJet(ilab);
412  jet2=jetUS->ClosestJet();
413 
414  if (!jet2) {
415  Printf("jet2 does not exist, returning");
416  continue;
417  }
418 
419  for(Int_t j=0; j<jetContPart->GetNJets(); j++) {
420 
421  jet3 = jetContPart->GetJet(j);
422  if(!jet3) continue;
423  if(jet3->GetLabel()==jet2->GetLabel()) {
424  jfound++;
425  if(jfound==1) jlab = j;
426  }
427  }
428  if(jlab==-1) continue;
429  jet3=jetContPart->GetJet(jlab);
430  if(!jet3){
431  Printf("jet3 does not exist, returning");
432  continue;
433  }
434  }
435  if(!(fJetShapeSub==kConstSub)) jet3 = jet1->ClosestJet();
436  if (!jet3) {
437  Printf("jet3 does not exist, returning");
438  continue;
439  }
440 
441 
442  fh2ResponseUW->Fill(jet1->Pt(),jet3->Pt());
443 
444 
445  }
446 
447 
448  if (fJetShapeType == kGenOnTheFly){
449  const AliEmcalPythiaInfo *partonsInfo = 0x0;
450  partonsInfo = GetPythiaInfo();
451  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
452  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
453  kWeight=partonsInfo->GetPythiaEventWeight();
454  fh2ResponseW->Fill(jet1->Pt(),jet1->Pt(),kWeight);
455 
456  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
457  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
458  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
459  if(dRp1 < fRMatching) {
460  fShapesVar[0] = partonsInfo->GetPartonFlag6();
461  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
462  }
463  else {
464  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
465  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
466  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
467  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
468  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
469  if(dRp1 < fRMatching) {
470  fShapesVar[0] = partonsInfo->GetPartonFlag7();
471  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
472  }
473  else fShapesVar[0]=0;
474  }
475  }
476 
477 
478 
479 
480  Double_t ptSubtracted = 0;
481  if (fJetShapeSub==kConstSub) ptSubtracted= jet1->Pt();
482 
483  else if (fJetShapeSub==kDerivSub) {
484  ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
485  }
486 
487  else if (fJetShapeSub==kNoSub) {
488  if ((fJetShapeType==kData) || (fJetShapeType==kDetEmbPartPythia)) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
489  else if ((fJetShapeType==kPythiaDef) || (fJetShapeType==kMCTrue) || (fJetShapeType==kGenOnTheFly)) ptSubtracted= jet1->Pt();
490  }
491 
492  //Printf("ptSubtracted=%f,fPtThreshold =%f ", ptSubtracted, fPtThreshold);
493  if (ptSubtracted < fPtThreshold) continue;
494 
495  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
496  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
497 
498 
499  fShapesVar[1] = ptSubtracted;
500  fShapesVar[2] = GetJetpTD(jet1,0);
501  fShapesVar[3] = GetJetMass(jet1,0);
502  fShapesVar[4] = GetJetAngularity(jet1,0);
503  fShapesVar[5] = GetJetCircularity(jet1,0);
504  fShapesVar[6] = 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  massMatch=GetJetMass(jet3,kMatched);
517  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
518  angulMatch=GetJetAngularity(jet3, kMatched);
519  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  massMatch=GetJetMass(jet3,kMatched);
530  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
531  angulMatch=GetJetAngularity(jet3, kMatched);
532  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 
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  fShapesVar[14] = ntracksEvt;
563  fShapesVar[15] = rhoVal;
564  fShapesVar[16] = rhoMassVal;
565  fShapesVar[17] = jet1->Pt();
566 
567 
568  fTreeObservableTagging->Fill();
569 
570 
571 
572 
573 
574 
575  }
576 
577  }
578 
579  return kTRUE;
580 }
581 
582 //________________________________________________________________________
584  //calc subtracted jet mass
585  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
587  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
588  else
589  return jet->M();
590 }
591 
592 //________________________________________________________________________
594 
595  AliJetContainer *jetCont = GetJetContainer(jetContNb);
596  if (!jet->GetNumberOfTracks())
597  return 0;
598  Double_t den=0.;
599  Double_t num = 0.;
600  AliVParticle *vp1 = 0x0;
601  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
602  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
603 
604  if (!vp1){
605  Printf("AliVParticle associated to constituent not found");
606  continue;
607  }
608 
609  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
610  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
611  Double_t dr = TMath::Sqrt(dr2);
612  num=num+vp1->Pt()*dr;
613  den=den+vp1->Pt();
614  }
615  return num/den;
616 }
617 
618 //________________________________________________________________________
620 
621  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
624  else
625  return Angularity(jet, jetContNb);
626 
627 }
628 
629 
630 //________________________________________________________________________
632 
633  AliJetContainer *jetCont = GetJetContainer(jetContNb);
634  if (!jet->GetNumberOfTracks())
635  return 0;
636  Double_t den=0.;
637  Double_t num = 0.;
638  AliVParticle *vp1 = 0x0;
639  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
640  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
641 
642  if (!vp1){
643  Printf("AliVParticle associated to constituent not found");
644  continue;
645  }
646 
647  num=num+vp1->Pt()*vp1->Pt();
648  den=den+vp1->Pt();
649  }
650  return TMath::Sqrt(num)/den;
651 }
652 
653 //________________________________________________________________________
655  //calc subtracted jet mass
656  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
658  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
659  else
660  return PTD(jet, jetContNb);
661 
662 }
663 
664 //_____________________________________________________________________________
666 
667  AliJetContainer *jetCont = GetJetContainer(jetContNb);
668  if (!jet->GetNumberOfTracks())
669  return 0;
670  Double_t mxx = 0.;
671  Double_t myy = 0.;
672  Double_t mxy = 0.;
673  int nc = 0;
674  Double_t sump2 = 0.;
675  Double_t pxjet=jet->Px();
676  Double_t pyjet=jet->Py();
677  Double_t pzjet=jet->Pz();
678 
679 
680  //2 general normalized vectors perpendicular to the jet
681  TVector3 ppJ1(pxjet, pyjet, pzjet);
682  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
683  ppJ3.SetMag(1.);
684  TVector3 ppJ2(-pyjet, pxjet, 0);
685  ppJ2.SetMag(1.);
686  AliVParticle *vp1 = 0x0;
687  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
688  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
689 
690  if (!vp1){
691  Printf("AliVParticle associated to constituent not found");
692  continue;
693  }
694 
695  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
696 
697  //local frame
698  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
699  TVector3 pPerp = pp - pLong;
700  //projection onto the two perpendicular vectors defined above
701 
702  Float_t ppjX = pPerp.Dot(ppJ2);
703  Float_t ppjY = pPerp.Dot(ppJ3);
704  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
705  if(ppjT<=0) return 0;
706 
707  mxx += (ppjX * ppjX / ppjT);
708  myy += (ppjY * ppjY / ppjT);
709  mxy += (ppjX * ppjY / ppjT);
710  nc++;
711  sump2 += ppjT;}
712 
713  if(nc<2) return 0;
714  if(sump2==0) return 0;
715  // Sphericity Matrix
716  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
717  TMatrixDSym m0(2,ele);
718 
719  // Find eigenvectors
720  TMatrixDSymEigen m(m0);
721  TVectorD eval(2);
722  TMatrixD evecm = m.GetEigenVectors();
723  eval = m.GetEigenValues();
724  // Largest eigenvector
725  int jev = 0;
726  // cout<<eval[0]<<" "<<eval[1]<<endl;
727  if (eval[0] < eval[1]) jev = 1;
728  TVectorD evec0(2);
729  // Principle axis
730  evec0 = TMatrixDColumn(evecm, jev);
731  Double_t compx=evec0[0];
732  Double_t compy=evec0[1];
733  TVector2 evec(compx, compy);
734  Double_t circ=0;
735  if(jev==1) circ=2*eval[0];
736  if(jev==0) circ=2*eval[1];
737 
738  return circ;
739 
740 
741 
742 }
743 
744 
745 
746 
747 //________________________________________________________________________
749  //calc subtracted jet mass
750 
751  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
754  else
755  return Circularity(jet, jetContNb);
756 
757 }
758 
759 //________________________________________________________________________
761 
762  AliJetContainer *jetCont = GetJetContainer(jetContNb);
763  if (!jet->GetNumberOfTracks())
764  return 0;
765  Double_t den=0.;
766  Double_t num = 0.;
767  AliVParticle *vp1 = 0x0;
768  AliVParticle *vp2 = 0x0;
769  std::vector<int> ordindex;
770  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
771  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
772  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
773 
774  if(ordindex.size()<2) return -1;
775 
776  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
777  if (!vp1){
778  Printf("AliVParticle associated to Leading constituent not found");
779  return -1;
780  }
781 
782  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
783  if (!vp2){
784  Printf("AliVParticle associated to Subleading constituent not found");
785  return -1;
786  }
787 
788 
789  num=vp1->Pt();
790  den=vp2->Pt();
791  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
792 
793 return num-den;
794 }
795 
796 //________________________________________________________________________
798  //calc subtracted jet mass
799 
800  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
803  else
804  return LeSub(jet, jetContNb);
805 
806 }
807 
808 //________________________________________________________________________
810  //calc subtracted jet mass
811 
812  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
815  else
816  return jet->GetNumberOfTracks();
817 
818  }
819 
820 
821 //______________________________________________________________________________
823 
824  AliJetContainer *jetCont = GetJetContainer(jetContNb);
825  if (!jet->GetNumberOfTracks())
826  return 0;
827  Double_t mxx = 0.;
828  Double_t myy = 0.;
829  Double_t mxy = 0.;
830  int nc = 0;
831  Double_t sump2 = 0.;
832 
833  AliVParticle *vp1 = 0x0;
834  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
835  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
836 
837  if (!vp1){
838  Printf("AliVParticle associated to constituent not found");
839  continue;
840  }
841 
842  Double_t ppt=vp1->Pt();
843  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
844 
845  Double_t deta = vp1->Eta()-jet->Eta();
846  mxx += ppt*ppt*deta*deta;
847  myy += ppt*ppt*dphi*dphi;
848  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
849  nc++;
850  sump2 += ppt*ppt;
851 
852  }
853  if(nc<2) return 0;
854  if(sump2==0) return 0;
855  // Sphericity Matrix
856  Double_t ele[4] = {mxx , mxy , mxy , myy };
857  TMatrixDSym m0(2,ele);
858 
859  // Find eigenvectors
860  TMatrixDSymEigen m(m0);
861  TVectorD eval(2);
862  TMatrixD evecm = m.GetEigenVectors();
863  eval = m.GetEigenValues();
864  // Largest eigenvector
865  int jev = 0;
866  // cout<<eval[0]<<" "<<eval[1]<<endl;
867  if (eval[0] < eval[1]) jev = 1;
868  TVectorD evec0(2);
869  // Principle axis
870  evec0 = TMatrixDColumn(evecm, jev);
871  Double_t compx=evec0[0];
872  Double_t compy=evec0[1];
873  TVector2 evec(compx, compy);
874  Double_t sig=0;
875  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
876  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
877 
878  return sig;
879 
880 }
881 
882 //________________________________________________________________________
884  //calc subtracted jet mass
885 
886  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
889  else
890  return Sigma2(jet, jetContNb);
891 
892 }
893 
894 //________________________________________________________________________
896 
898  TClonesArray *tracksArray = partCont->GetArray();
899 
900  if(!partCont || !tracksArray) return -99999;
901  AliAODTrack *track = 0x0;
902  AliEmcalParticle *emcPart = 0x0;
903 
904 
905  TList *trackList = new TList();
906  Int_t triggers[100];
907  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
908  Int_t iTT = 0;
909 
910  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
911 
912 
913  if (fJetShapeSub == kConstSub){
914  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
915  if (!emcPart) continue;
916  if(TMath::Abs(emcPart->Eta())>0.9) continue;
917  if (emcPart->Pt()<0.15) continue;
918 
919  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
920  trackList->Add(emcPart);
921  triggers[iTT] = iTrack;
922  iTT++;
923  }
924  }
925  else{
926  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
927  if (!track) continue;
928  if(TMath::Abs(track->Eta())>0.9) continue;
929  if (track->Pt()<0.15) continue;
930  if (!(track->TestFilterBit(768))) continue;
931 
932  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
933  trackList->Add(track);
934  triggers[iTT] = iTrack;
935  iTT++;
936 
937  }
938  }
939  }
940 
941  if (iTT == 0) return -99999;
942  Int_t nbRn = 0, index = 0 ;
943  TRandom3* random = new TRandom3(0);
944  nbRn = random->Integer(iTT);
945 
946  index = triggers[nbRn];
947  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
948  return index;
949 
950 }
951 
952 //__________________________________________________________________________________
954 
955  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
956  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
957  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
958  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
959  double dphi = mphi-vphi;
960  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
961  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
962  return dphi;//dphi in [-Pi, Pi]
963 }
964 
965 
966 //________________________________________________________________________
968  //
969  // retrieve event objects
970  //
972  return kFALSE;
973 
974  return kTRUE;
975 }
976 
977 //_______________________________________________________________________
979 {
980  // Called once at the end of the analysis.
981 
982  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
983  // if (!fTreeObservableTagging){
984  // Printf("ERROR: fTreeObservableTagging not available");
985  // return;
986  // }
987 
988 }
989 
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Double_t GetFirstOrderSubtractedAngularity() const
Double_t Area() const
Definition: AliEmcalJet.h:123
const TString & GetRhoName() const
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:219
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPythiaEventWeight() const
Int_t GetLabel() const
Definition: AliEmcalJet.h:117
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:132
TString kData
Declare data MC or deltaAOD.
Double_t Px() const
Definition: AliEmcalJet.h:99
AliParticleContainer * GetParticleContainer(Int_t i=0) const
const TString & GetRhoMassName() 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:102
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:153
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:101
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:253
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPartonPt7() const
Double_t M() const
Definition: AliEmcalJet.h:113
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