AliPhysics  4646b6b (4646b6b)
 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 //just testing
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"
45 #include "AliAnalysisTaskEmcalQGTagging.h"
46 
47 using std::cout;
48 using std::endl;
49 
50 ClassImp(AliAnalysisTaskEmcalQGTagging)
51 
52 //________________________________________________________________________
53 AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging() :
54  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalQGTagging", kTRUE),
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  fHLundIterative(0x0),
86  fNbOfConstvspT(0x0),
87  fTreeObservableTagging(0)
88 
89 {
90  for(Int_t i=0;i<17;i++){
91  fShapesVar[i]=0;}
92  SetMakeGeneralHistograms(kTRUE);
93  DefineOutput(1, TList::Class());
94  DefineOutput(2, TTree::Class());
95 }
96 
97 //________________________________________________________________________
98 AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) :
99  AliAnalysisTaskEmcalJet(name, kTRUE),
100  fContainer(0),
101  fMinFractionShared(0),
102  fJetShapeType(kData),
103  fJetShapeSub(kNoSub),
104  fJetSelection(kInclusive),
105  fPtThreshold(-9999.),
106  fRMatching(0.2),
107  fSelectedShapes(0),
108  fminpTTrig(20.),
109  fmaxpTTrig(50.),
110  fangWindowRecoil(0.6),
111  fSemigoodCorrect(0),
112  fHolePos(0),
113  fHoleWidth(0),
114  fCentSelectOn(kTRUE),
115  fCentMin(0),
116  fCentMax(10),
117  fOneConstSelectOn(kFALSE),
118  fDerivSubtrOrder(0),
119  fh2ResponseUW(0x0),
120  fh2ResponseW(0x0),
121  fPhiJetCorr6(0x0),
122  fPhiJetCorr7(0x0),
123  fEtaJetCorr6(0x0),
124  fEtaJetCorr7(0x0),
125  fPtJetCorr(0x0),
126  fPtJet(0x0),
127  fhpTjetpT(0x0),
128  fhPt(0x0),
129  fhPhi(0x0),
130  fHLundIterative(0x0),
131  fNbOfConstvspT(0x0),
132  fTreeObservableTagging(0)
133 
134 {
135  // Standard constructor.
136  for(Int_t i=0;i<17;i++){
137  fShapesVar[i]=0;}
138  SetMakeGeneralHistograms(kTRUE);
139 
140  DefineOutput(1, TList::Class());
141  DefineOutput(2, TTree::Class());
142 }
143 
144 //________________________________________________________________________
145 AliAnalysisTaskEmcalQGTagging::~AliAnalysisTaskEmcalQGTagging()
146 {
147  // Destructor.
148 }
149 
150 //________________________________________________________________________
151  void AliAnalysisTaskEmcalQGTagging::UserCreateOutputObjects()
152 {
153  // Create user output.
154 
155  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
156 
157  Bool_t oldStatus = TH1::AddDirectoryStatus();
158  TH1::AddDirectory(kFALSE);
159 
160 
161 
162 
163  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
164  fOutput->Add(fh2ResponseUW);
165  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
166  fOutput->Add(fh2ResponseW);
167  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
168  fOutput->Add(fPhiJetCorr6);
169  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
170  fOutput->Add(fEtaJetCorr6);
171 
172  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
173  fOutput->Add(fPhiJetCorr7);
174  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
175  fOutput->Add(fEtaJetCorr7);
176 
177  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
178  fOutput->Add(fPtJetCorr);
179  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
180  fOutput->Add(fPtJet);
181 
182  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
183  fOutput->Add(fhpTjetpT);
184  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
185  fOutput->Add(fhPt);
186  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
187  fOutput->Add(fhPhi);
188 
189  //log(1/theta),log(z*theta),jetpT,algo//
190  const Int_t dimSpec = 5;
191  const Int_t nBinsSpec[5] = {50,50,10,3,10};
192  const Double_t lowBinSpec[5] = {0.0,-10, 0,0,0};
193  const Double_t hiBinSpec[5] = {5.0, 0,200,3,10};
194  fHLundIterative = new THnSparseF("fHLundIterative",
195  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo]",
196  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
197  fOutput->Add(fHLundIterative);
198 
199  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
200  fOutput->Add(fNbOfConstvspT);
201 
202  // =========== Switch on Sumw2 for all histos ===========
203  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
204  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
205  if (h1){
206  h1->Sumw2();
207  continue;
208  }
209  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
210  if(hn)hn->Sumw2();
211  }
212 
213 
214  TH1::AddDirectory(oldStatus);
215  const Int_t nVar = 12;
216  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
217  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
218 
219 
220  TString *fShapesVarNames = new TString [nVar];
221 
222  fShapesVarNames[0] = "partonCode";
223  fShapesVarNames[1] = "ptJet";
224  fShapesVarNames[2] = "ptDJet";
225  fShapesVarNames[3] = "phiJet";
226  // fShapesVarNames[4] = "nbOfConst";
227  fShapesVarNames[4] = "angularity";
228  //fShapesVarNames[5] = "circularity";
229  fShapesVarNames[5] = "lesub";
230  //fShapesVarNames[7] = "coronna";
231 
232  fShapesVarNames[6] = "ptJetMatch";
233  fShapesVarNames[7] = "ptDJetMatch";
234  fShapesVarNames[8] = "phiJetMatch";
235  // fShapesVarNames[12] = "nbOfConstMatch";
236  fShapesVarNames[9] = "angularityMatch";
237  //fShapesVarNames[12] = "circularityMatch";
238  fShapesVarNames[10] = "lesubMatch";
239  //fShapesVarNames[14] = "coronnaMatch";
240  fShapesVarNames[11]="weightPythia";
241  //fShapesVarNames[14]="ntrksEvt";
242  //fShapesVarNames[16]="rhoVal";
243  //fShapesVarNames[17]="rhoMassVal";
244  //fShapesVarNames[12]="ptUnsub";
245 
246  for(Int_t ivar=0; ivar < nVar; ivar++){
247  cout<<"looping over variables"<<endl;
248  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
249 
250  PostData(1,fOutput);
251  PostData(2,fTreeObservableTagging);
252 
253  delete [] fShapesVarNames;
254 }
255 
256 //________________________________________________________________________
257 Bool_t AliAnalysisTaskEmcalQGTagging::Run()
258 {
259  // Run analysis code here, if needed. It will be executed before FillHistograms().
260 
261  return kTRUE;
262 }
263 
264 //________________________________________________________________________
265 Bool_t AliAnalysisTaskEmcalQGTagging::FillHistograms()
266 {
267  // Fill histograms.
268  //cout<<"base container"<<endl;
269  AliEmcalJet* jet1 = NULL;
270  AliJetContainer *jetCont = GetJetContainer(0);
271 
272  Float_t kWeight=1;
273  if (fCentSelectOn)
274  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
275 
276  AliAODTrack *triggerHadron = 0x0;
277 
278  if (fJetSelection == kRecoil) {
279  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
280  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
281 
282 
283  if (triggerHadronLabel==-99999) {
284  //Printf ("Trigger Hadron not found, return");
285  return 0;}
286 
287  AliTrackContainer *PartCont =NULL;
288  AliParticleContainer *PartContMC=NULL;
289 
290  if (fJetShapeSub==kConstSub){
291  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(1);
292  else PartCont = GetTrackContainer(1);
293  }
294  else{
295  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(0);
296  else PartCont = GetTrackContainer(0);
297  }
298  TClonesArray *TrackArray = NULL;
299  TClonesArray *TrackArrayMC = NULL;
300  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
301  else TrackArray = PartCont->GetArray();
302  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) triggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(triggerHadronLabel));
303  else triggerHadron = static_cast<AliAODTrack*>(TrackArray->At(triggerHadronLabel));
304 
305 
306 
307 
308 
309  if (!triggerHadron) {
310  //Printf("No Trigger hadron with the found label!!");
311  return 0;
312  }
313 
314  if(fSemigoodCorrect){
315  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
316  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
317  return 0;}
318  }
319 
320  fhPt->Fill(triggerHadron->Pt());
321 
322  }
323 
324 
325  AliParticleContainer *partContAn = GetParticleContainer(0);
326  TClonesArray *trackArrayAn = partContAn->GetArray();
327  Int_t ntracksEvt = trackArrayAn->GetEntriesFast();
328 
329  Float_t rhoVal=0, rhoMassVal = 0.;
330  if(jetCont) {
331 
332  jetCont->ResetCurrentID();
333  if ((fJetShapeSub==kConstSub) || (fJetShapeSub==kDerivSub)){
334  //rho
335  AliRhoParameter* rhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoSparseR020"));
336  if (!rhoParam) {
337  Printf("%s: Could not retrieve rho %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoName().Data());
338  } else rhoVal = rhoParam->GetVal();
339  //rhom
340  AliRhoParameter* rhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoMassSparseR020"));
341  if (!rhomParam) {
342  Printf("%s: Could not retrieve rho_m %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoMassName().Data());
343  } else rhoMassVal = rhomParam->GetVal();
344  }
345 
346  while((jet1 = jetCont->GetNextAcceptJet())) {
347  if (!jet1) continue;
348  AliEmcalJet* jet2 = 0x0;
349  AliEmcalJet* jet3 = 0x0;
350  fPtJet->Fill(jet1->Pt());
351  AliEmcalJet *jetUS = NULL;
352  Int_t ifound=0, jfound=0;
353  Int_t ilab=-1, jlab=-1;
354 
355  if(fSemigoodCorrect && (fJetSelection != kRecoil)){
356  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
357  if(TMath::Abs(disthole)<fHoleWidth){
358  continue;}
359  }
360 
361  Float_t dphiRecoil = 0.;
362  if (fJetSelection == kRecoil){
363  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
364  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
365  // Printf("Recoil jets back to back not found! continuing");
366  continue;
367  }
368 
369  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
370  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
371  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
372 
373  }
374 
375 
376  fShapesVar[0] = 0.;
377  if(fJetShapeType == kDetEmbPartPythia){
378  AliJetContainer *jetContTrue = GetJetContainer(1);
379  AliJetContainer *jetContUS = GetJetContainer(2);
380 
381  if(fJetShapeSub==kConstSub){
382  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
383  jetUS = jetContUS->GetJet(i);
384  if(jetUS->GetLabel()==jet1->GetLabel()) {
385  ifound++;
386  if(ifound==1) ilab = i;
387  }
388  }
389  if(ilab==-1) continue;
390  jetUS=jetContUS->GetJet(ilab);
391  jet2=jetUS->ClosestJet();
392  }
393 
394  if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet();
395  if (!jet2) {
396  Printf("jet2 does not exist, returning");
397  continue;
398  }
399 
400  AliJetContainer *jetContPart=GetJetContainer(3);
401  jet3=jet2->ClosestJet();
402 
403  if(!jet3){
404  Printf("jet3 does not exist, returning");
405  continue;
406  }
407  cout<<"jet 3 exists"<<jet3->Pt()<<endl;
408 
409 
410  fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt());
411 
412  Double_t fraction=0;
413  if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1);
414  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
415  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
416 
417  if(fraction<fMinFractionShared) continue;
418  //InputEvent()->Print();
419 
420  }
421 
422 
423 
424  if (fJetShapeType == kPythiaDef){
425 
426  AliJetContainer *jetContTrue = GetJetContainer(1);
427  AliJetContainer *jetContUS = GetJetContainer(2);
428  AliJetContainer *jetContPart = GetJetContainer(3);
429 
430  if(fJetShapeSub==kConstSub){
431 
432  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
433  jetUS = jetContUS->GetJet(i);
434  if(jetUS->GetLabel()==jet1->GetLabel()) {
435  ifound++;
436  if(ifound==1) ilab = i;
437  }
438  }
439  if(ilab==-1) continue;
440  jetUS=jetContUS->GetJet(ilab);
441  jet2=jetUS->ClosestJet();
442 
443  if (!jet2) {
444  Printf("jet2 does not exist, returning");
445  continue;
446  }
447 
448  for(Int_t j=0; j<jetContPart->GetNJets(); j++) {
449 
450  jet3 = jetContPart->GetJet(j);
451  if(!jet3) continue;
452  if(jet3->GetLabel()==jet2->GetLabel()) {
453  jfound++;
454  if(jfound==1) jlab = j;
455  }
456  }
457  if(jlab==-1) continue;
458  jet3=jetContPart->GetJet(jlab);
459  if(!jet3){
460  Printf("jet3 does not exist, returning");
461  continue;
462  }
463  }
464  if(!(fJetShapeSub==kConstSub)) jet3 = jet1->ClosestJet();
465  if (!jet3) {
466  Printf("jet3 does not exist, returning");
467  continue;
468  }
469 
470 
471  fh2ResponseUW->Fill(jet1->Pt(),jet3->Pt());
472 
473 
474  }
475 
476 
477  if (fJetShapeType == kGenOnTheFly){
478  const AliEmcalPythiaInfo *partonsInfo = 0x0;
479  partonsInfo = GetPythiaInfo();
480  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
481  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
482  kWeight=partonsInfo->GetPythiaEventWeight();
483  fh2ResponseW->Fill(jet1->Pt(),jet1->Pt(),kWeight);
484 
485  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
486  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
487  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
488  if(dRp1 < fRMatching) {
489  fShapesVar[0] = partonsInfo->GetPartonFlag6();
490  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
491  }
492  else {
493  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
494  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
495  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
496  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
497  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
498  if(dRp1 < fRMatching) {
499  fShapesVar[0] = partonsInfo->GetPartonFlag7();
500  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
501  }
502  else fShapesVar[0]=0;
503  }
504  }
505 
506 
507 
508 
509  Double_t ptSubtracted = 0;
510  if (fJetShapeSub==kConstSub) ptSubtracted= jet1->Pt();
511 
512  else if (fJetShapeSub==kDerivSub) {
513  ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
514  }
515 
516  else if (fJetShapeSub==kNoSub) {
517  if ((fJetShapeType==kData) || (fJetShapeType==kDetEmbPartPythia)) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
518  else if ((fJetShapeType==kPythiaDef) || (fJetShapeType==kMCTrue) || (fJetShapeType==kGenOnTheFly)) ptSubtracted= jet1->Pt();
519  }
520 
521  //Printf("ptSubtracted=%f,fPtThreshold =%f ", ptSubtracted, fPtThreshold);
522  if (ptSubtracted < fPtThreshold) continue;
523 
524  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
525  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
526 
527 
528  fShapesVar[1] = ptSubtracted;
529  fShapesVar[2] = GetJetpTD(jet1,0);
530  fShapesVar[3] =jet1->Phi();
531  if(fJetShapeType==kData) fShapesVar[3]=RelativePhi(triggerHadron->Phi(), jet1->Phi());
532  //GetJetMass(jet1,0);
533  fShapesVar[4] = GetJetAngularity(jet1,0);
534  //fShapesVar[5] = GetJetCircularity(jet1,0);
535  fShapesVar[5] = GetJetLeSub(jet1,0);
536  //fShapesVar[6] = GetJetCoronna(jet1,0);
537  RecursiveParents(jet1,jetCont,0);
538  RecursiveParents(jet1,jetCont,1);
539  RecursiveParents(jet1,jetCont,2);
540 
541  Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0., coronnaMatch=0;
542  Int_t kMatched = 0;
543 
544  if (fJetShapeType==kPythiaDef) {
545  kMatched =1;
546  if(fJetShapeSub==kConstSub) kMatched = 3;
547 
548  ptMatch=jet3->Pt();
549  ptDMatch=GetJetpTD(jet3, kMatched);
550  massMatch=jet3->Phi();
551  // GetJetMass(jet3,kMatched);
552  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
553  angulMatch=GetJetAngularity(jet3, kMatched);
554  //circMatch=GetJetCircularity(jet3, kMatched);
555  lesubMatch=GetJetLeSub(jet3, kMatched);
556  //coronnaMatch=GetJetCoronna(jet3,kMatched);
557  //sigma2Match = GetSigma2(jet2, kMatched);
558  }
559 
560  if (fJetShapeType==kDetEmbPartPythia) {
561  if(fJetShapeSub==kConstSub) kMatched = 3;
562  if(fJetShapeSub==kDerivSub) kMatched = 2;
563  ptMatch=jet3->Pt();
564  ptDMatch=GetJetpTD(jet3, kMatched);
565  massMatch=jet3->Phi();
566  //GetJetMass(jet3,kMatched);
567  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
568  angulMatch=GetJetAngularity(jet3, kMatched);
569  // circMatch=GetJetCircularity(jet3, kMatched);
570  lesubMatch=GetJetLeSub(jet3, kMatched);
571  //coronnaMatch = GetJetCoronna(jet3, kMatched);
572 
573  }
574 
575 
576 
577  if (fJetShapeType == kMCTrue || fJetShapeType == kData || fJetShapeType == kGenOnTheFly) {
578  kMatched = 0;
579  ptMatch=0.;
580  ptDMatch=0.;
581  massMatch=0.;
582  //constMatch=0.;
583  angulMatch=0.;
584  // circMatch=0.;
585  lesubMatch=0.;
586  //coronnaMatch =0.;
587 
588  }
589 
590 
591 
592  fShapesVar[6] = ptMatch;
593  fShapesVar[7] = ptDMatch;
594  fShapesVar[8] = massMatch;
595  fShapesVar[9] = angulMatch;
596  //fShapesVar[12] = circMatch;
597  fShapesVar[10] = lesubMatch;
598  // fShapesVar[14] = coronnaMatch;
599  fShapesVar[11] = kWeight;
600  //fShapesVar[16] = ntracksEvt;
601  // fShapesVar[16] = rhoVal;
602  //fShapesVar[17] = rhoMassVal;
603  //fShapesVar[16] = jet1->Pt();
604 
605 
606  fTreeObservableTagging->Fill();
607 
608 
609 
610 
611 
612 
613  }
614 
615  }
616 
617  return kTRUE;
618 }
619 
620 //________________________________________________________________________
621 Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0){
622  //calc subtracted jet mass
623  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
624  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtracted();
625  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
626  else
627  return jet->M();
628 }
629 
630 //________________________________________________________________________
631 Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
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  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
648  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
649  Double_t dr = TMath::Sqrt(dr2);
650  num=num+vp1->Pt()*dr;
651  den=den+vp1->Pt();
652  }
653  return num/den;
654 }
655 
656 //________________________________________________________________________
657 Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0){
658 
659  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
660  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedAngularity();
661  else return jet->GetShapeProperties()->GetSecondOrderSubtractedAngularity();
662  else
663  return Angularity(jet, jetContNb);
664 
665 }
666 
667 //____________________________________________________________________________
668 
669 Float_t AliAnalysisTaskEmcalQGTagging::Coronna(AliEmcalJet *jet, Int_t jetContNb = 0){
670 
671  AliTrackContainer *PartCont = NULL;
672  AliParticleContainer *PartContMC = NULL;
673 
674 
675  if (fJetShapeSub==kConstSub){
676  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(1);
677  else PartCont = GetTrackContainer(1);
678  }
679  else{
680  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(0);
681  else PartCont = GetTrackContainer(0);
682  }
683 
684  TClonesArray *TracksArray = NULL;
685  TClonesArray *TracksArrayMC = NULL;
686 
687  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
688  else TracksArray = PartCont->GetArray();
689 
690  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
691  if(!PartContMC || !TracksArrayMC) return -2;
692  }
693  else {
694  if(!PartCont || !TracksArray) return -2;
695  }
696 
697 
698  AliAODTrack *Track = 0x0;
699  Float_t sumpt=0;
700  Int_t NTracks=0;
701  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
702  else NTracks = TracksArray->GetEntriesFast();
703 
704  for(Int_t i=0; i < NTracks; i++){
705  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
706  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
707  if (!Track) continue;
708  if(TMath::Abs(Track->Eta())>0.9) continue;
709  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
710  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
711  Double_t dr = TMath::Sqrt(dr2);
712  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
713  }
714  }
715  else{
716  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
717  if (!Track) continue;
718  if(TMath::Abs(Track->Eta())>0.9) continue;
719  if (Track->Pt()<0.15) continue;
720  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
721  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
722  Double_t dr = TMath::Sqrt(dr2);
723  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
724 
725  }
726  }
727  }
728 
729 
730 
731  return sumpt;
732 
733 
734 
735 
736 }
737 
738 //________________________________________________________________________
739 Float_t AliAnalysisTaskEmcalQGTagging::GetJetCoronna(AliEmcalJet *jet, Int_t jetContNb = 0){
740 
741  if((fJetShapeSub==kDerivSub) && (jetContNb==0)) return -2;
742  else
743  return Coronna(jet, jetContNb);
744 
745 }
746 
747 
748 
749 
750 
751 //________________________________________________________________________
752 Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
753 
754  AliJetContainer *jetCont = GetJetContainer(jetContNb);
755  if (!jet->GetNumberOfTracks())
756  return 0;
757  Double_t den=0.;
758  Double_t num = 0.;
759  AliVParticle *vp1 = 0x0;
760  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
761  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
762 
763  if (!vp1){
764  Printf("AliVParticle associated to constituent not found");
765  continue;
766  }
767 
768  num=num+vp1->Pt()*vp1->Pt();
769  den=den+vp1->Pt();
770  }
771  return TMath::Sqrt(num)/den;
772 }
773 
774 //________________________________________________________________________
775 Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0){
776  //calc subtracted jet mass
777  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
778  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedpTD();
779  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
780  else
781  return PTD(jet, jetContNb);
782 
783 }
784 
785 //_____________________________________________________________________________
786 Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
787 
788  AliJetContainer *jetCont = GetJetContainer(jetContNb);
789  if (!jet->GetNumberOfTracks())
790  return 0;
791  Double_t mxx = 0.;
792  Double_t myy = 0.;
793  Double_t mxy = 0.;
794  int nc = 0;
795  Double_t sump2 = 0.;
796  Double_t pxjet=jet->Px();
797  Double_t pyjet=jet->Py();
798  Double_t pzjet=jet->Pz();
799 
800 
801  //2 general normalized vectors perpendicular to the jet
802  TVector3 ppJ1(pxjet, pyjet, pzjet);
803  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
804  ppJ3.SetMag(1.);
805  TVector3 ppJ2(-pyjet, pxjet, 0);
806  ppJ2.SetMag(1.);
807  AliVParticle *vp1 = 0x0;
808  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
809  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
810 
811  if (!vp1){
812  Printf("AliVParticle associated to constituent not found");
813  continue;
814  }
815 
816  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
817 
818  //local frame
819  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
820  TVector3 pPerp = pp - pLong;
821  //projection onto the two perpendicular vectors defined above
822 
823  Float_t ppjX = pPerp.Dot(ppJ2);
824  Float_t ppjY = pPerp.Dot(ppJ3);
825  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
826  if(ppjT<=0) return 0;
827 
828  mxx += (ppjX * ppjX / ppjT);
829  myy += (ppjY * ppjY / ppjT);
830  mxy += (ppjX * ppjY / ppjT);
831  nc++;
832  sump2 += ppjT;}
833 
834  if(nc<2) return 0;
835  if(sump2==0) return 0;
836  // Sphericity Matrix
837  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
838  TMatrixDSym m0(2,ele);
839 
840  // Find eigenvectors
841  TMatrixDSymEigen m(m0);
842  TVectorD eval(2);
843  TMatrixD evecm = m.GetEigenVectors();
844  eval = m.GetEigenValues();
845  // Largest eigenvector
846  int jev = 0;
847  // cout<<eval[0]<<" "<<eval[1]<<endl;
848  if (eval[0] < eval[1]) jev = 1;
849  TVectorD evec0(2);
850  // Principle axis
851  evec0 = TMatrixDColumn(evecm, jev);
852  Double_t compx=evec0[0];
853  Double_t compy=evec0[1];
854  TVector2 evec(compx, compy);
855  Double_t circ=0;
856  if(jev==1) circ=2*eval[0];
857  if(jev==0) circ=2*eval[1];
858 
859  return circ;
860 
861 
862 
863 }
864 
865 
866 
867 
868 //________________________________________________________________________
869 Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ){
870  //calc subtracted jet mass
871 
872  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
873  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedCircularity();
874  else return jet->GetShapeProperties()->GetSecondOrderSubtractedCircularity();
875  else
876  return Circularity(jet, jetContNb);
877 
878 }
879 
880 //________________________________________________________________________
881 Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
882 
883  AliJetContainer *jetCont = GetJetContainer(jetContNb);
884  if (!jet->GetNumberOfTracks())
885  return 0;
886  Double_t den=0.;
887  Double_t num = 0.;
888  AliVParticle *vp1 = 0x0;
889  AliVParticle *vp2 = 0x0;
890  std::vector<int> ordindex;
891  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
892  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
893  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
894 
895  if(ordindex.size()<2) return -1;
896 
897  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
898  if (!vp1){
899  Printf("AliVParticle associated to Leading constituent not found");
900  return -1;
901  }
902 
903  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
904  if (!vp2){
905  Printf("AliVParticle associated to Subleading constituent not found");
906  return -1;
907  }
908 
909 
910  num=vp1->Pt();
911  den=vp2->Pt();
912  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
913 
914 return num-den;
915 }
916 
917 //________________________________________________________________________
918 Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
919  //calc subtracted jet mass
920 
921  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
922  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedLeSub();
923  else return jet->GetShapeProperties()->GetSecondOrderSubtractedLeSub();
924  else
925  return LeSub(jet, jetContNb);
926 
927 }
928 
929 //________________________________________________________________________
930  Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0){
931  //calc subtracted jet mass
932 
933  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
934  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedConstituent();
935  else return jet->GetShapeProperties()->GetSecondOrderSubtractedConstituent();
936  else
937  return jet->GetNumberOfTracks();
938 
939  }
940 
941 
942 //______________________________________________________________________________
943 Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
944 
945  AliJetContainer *jetCont = GetJetContainer(jetContNb);
946  if (!jet->GetNumberOfTracks())
947  return 0;
948  Double_t mxx = 0.;
949  Double_t myy = 0.;
950  Double_t mxy = 0.;
951  int nc = 0;
952  Double_t sump2 = 0.;
953 
954  AliVParticle *vp1 = 0x0;
955  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
956  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
957 
958  if (!vp1){
959  Printf("AliVParticle associated to constituent not found");
960  continue;
961  }
962 
963  Double_t ppt=vp1->Pt();
964  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
965 
966  Double_t deta = vp1->Eta()-jet->Eta();
967  mxx += ppt*ppt*deta*deta;
968  myy += ppt*ppt*dphi*dphi;
969  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
970  nc++;
971  sump2 += ppt*ppt;
972 
973  }
974  if(nc<2) return 0;
975  if(sump2==0) return 0;
976  // Sphericity Matrix
977  Double_t ele[4] = {mxx , mxy , mxy , myy };
978  TMatrixDSym m0(2,ele);
979 
980  // Find eigenvectors
981  TMatrixDSymEigen m(m0);
982  TVectorD eval(2);
983  TMatrixD evecm = m.GetEigenVectors();
984  eval = m.GetEigenValues();
985  // Largest eigenvector
986  int jev = 0;
987  // cout<<eval[0]<<" "<<eval[1]<<endl;
988  if (eval[0] < eval[1]) jev = 1;
989  TVectorD evec0(2);
990  // Principle axis
991  evec0 = TMatrixDColumn(evecm, jev);
992  Double_t compx=evec0[0];
993  Double_t compy=evec0[1];
994  TVector2 evec(compx, compy);
995  Double_t sig=0;
996  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
997  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
998 
999  return sig;
1000 
1001 }
1002 
1003 //________________________________________________________________________
1004 Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0){
1005  //calc subtracted jet mass
1006 
1007  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1008  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedSigma2();
1009  else return jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2();
1010  else
1011  return Sigma2(jet, jetContNb);
1012 
1013 }
1014 
1015 //________________________________________________________________________
1016 Int_t AliAnalysisTaskEmcalQGTagging::SelectTrigger(Float_t minpT, Float_t maxpT){
1017 
1018  AliTrackContainer *PartCont = NULL;
1019  AliParticleContainer *PartContMC = NULL;
1020 
1021 
1022  if (fJetShapeSub==kConstSub){
1023  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(1);
1024  else PartCont = GetTrackContainer(1);
1025  }
1026  else{
1027  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(0);
1028  else PartCont = GetTrackContainer(0);
1029  }
1030 
1031  TClonesArray *TracksArray = NULL;
1032  TClonesArray *TracksArrayMC = NULL;
1033 
1034  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1035  else TracksArray = PartCont->GetArray();
1036 
1037  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
1038  if(!PartContMC || !TracksArrayMC) return -99999;
1039  }
1040  else {
1041  if(!PartCont || !TracksArray) return -99999;
1042  }
1043 
1044 
1045  AliAODTrack *Track = 0x0;
1046 
1047 
1048 
1049  TList *trackList = new TList();
1050  Int_t triggers[100];
1051  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1052  Int_t iTT = 0;
1053  Int_t NTracks=0;
1054  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1055  else NTracks = TracksArray->GetEntriesFast();
1056 
1057  for(Int_t i=0; i < NTracks; i++){
1058  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
1059  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1060  if (!Track) continue;
1061  if(TMath::Abs(Track->Eta())>0.9) continue;
1062  if (Track->Pt()<0.15) continue;
1063  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1064  triggers[iTT] = i;
1065  iTT++;
1066  }
1067  }
1068  }
1069  else{
1070  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1071  if (!Track) continue;
1072  if(TMath::Abs(Track->Eta())>0.9) continue;
1073  if (Track->Pt()<0.15) continue;
1074  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1075  triggers[iTT] = i;
1076  iTT++;
1077  }
1078  }
1079  }
1080  }
1081 
1082 
1083  if (iTT == 0) return -99999;
1084  Int_t nbRn = 0, index = 0 ;
1085  TRandom3* random = new TRandom3(0);
1086  nbRn = random->Integer(iTT);
1087  index = triggers[nbRn];
1088  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1089  return index;
1090 
1091 }
1092 
1093 //__________________________________________________________________________________
1094 Double_t AliAnalysisTaskEmcalQGTagging::RelativePhi(Double_t mphi,Double_t vphi){
1095 
1096  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1097  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1098  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1099  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1100  double dphi = mphi-vphi;
1101  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1102  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1103  return dphi;//dphi in [-Pi, Pi]
1104 }
1105 
1106 
1107 //_________________________________________________________________________
1108 void AliAnalysisTaskEmcalQGTagging::RecursiveParents(AliEmcalJet *fJet,AliJetContainer *fJetCont, Int_t ReclusterAlgo){
1109 
1110  std::vector<fastjet::PseudoJet> fInputVectors;
1111  fInputVectors.clear();
1112  fastjet::PseudoJet PseudoTracks;
1113  double xflagalgo=0;
1114  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1115 
1116  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1117  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1118  if (!fTrk) continue;
1119  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1120  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1121  fInputVectors.push_back(PseudoTracks);
1122 
1123  }
1124  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1125 
1126  if(ReclusterAlgo==0){ xflagalgo=0.5;
1127  jetalgo=fastjet::kt_algorithm ;}
1128 
1129  if(ReclusterAlgo==1){ xflagalgo=1.5;
1130  jetalgo=fastjet::cambridge_algorithm;}
1131  if(ReclusterAlgo==2){ xflagalgo=2.5;
1132  jetalgo=fastjet::antikt_algorithm;}
1133 
1134  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1135 
1136  try {
1137  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1138  std::vector<fastjet::PseudoJet> fOutputJets;
1139  fOutputJets.clear();
1140  fOutputJets=fClustSeqSA.inclusive_jets(0);
1141 
1142  fastjet::PseudoJet jj;
1143  fastjet::PseudoJet j1;
1144  fastjet::PseudoJet j2;
1145  jj=fOutputJets[0];
1146  double ndepth=0;
1147  while(jj.has_parents(j1,j2)){
1148  ndepth=ndepth+1;
1149  if(j1.perp() < j2.perp()) swap(j1,j2);
1150  double delta_R=j1.delta_R(j2);
1151  double z=j2.perp()/(j1.perp()+j2.perp());
1152  double y =log(1.0/delta_R);
1153  double lnpt_rel=log(z*delta_R);
1154  Double_t LundEntries[5] = {y,lnpt_rel,fOutputJets[0].perp(),xflagalgo,ndepth};
1155  fHLundIterative->Fill(LundEntries);
1156  jj=j1;}
1157 
1158 
1159 
1160 
1161  } catch (fastjet::Error) {
1162  AliError(" [w] FJ Exception caught.");
1163  //return -1;
1164  }
1165 
1166 
1167 
1168 
1169  return;
1170 
1171 
1172 }
1173 
1174 
1175 
1176 
1177 
1178 
1179 
1180 //________________________________________________________________________
1181 Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
1182  //
1183  // retrieve event objects
1184  //
1185  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
1186  return kFALSE;
1187 
1188  return kTRUE;
1189 }
1190 
1191 //_______________________________________________________________________
1192 void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *)
1193 {
1194  // Called once at the end of the analysis.
1195 
1196  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
1197  // if (!fTreeObservableTagging){
1198  // Printf("ERROR: fTreeObservableTagging not available");
1199  // return;
1200  // }
1201 
1202 }
1203 
AliAnalysisTaskEmcalQGTagging::SelectTrigger
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1016
AliAnalysisTaskEmcalQGTagging::Coronna
Float_t Coronna(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:669
AliEmcalJetShapeProperties::GetFirstOrderSubtractedAngularity
Double_t GetFirstOrderSubtractedAngularity() const
Definition: AliEmcalJetShapeProperties.h:62
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:130
AliJetContainer::GetRhoName
const TString & GetRhoName() const
Definition: AliJetContainer.h:163
Double_t
double Double_t
Definition: External.C:58
AliEmcalJetShapeProperties::GetSecondOrderSubtractedSigma2
Double_t GetSecondOrderSubtractedSigma2() const
Definition: AliEmcalJetShapeProperties.h:93
TH2F
Definition: External.C:236
AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1181
AliEmcalPythiaInfo::GetPartonEta6
Float_t GetPartonEta6() const
Definition: AliEmcalPythiaInfo.h:45
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:377
AliAnalysisTaskEmcalQGTagging::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskEmcalQGTagging.cxx:265
AliEmcalPythiaInfo::GetPartonEta7
Float_t GetPartonEta7() const
Definition: AliEmcalPythiaInfo.h:50
AliEmcalJetShapeProperties::GetSecondOrderSubtractedConstituent
Double_t GetSecondOrderSubtractedConstituent() const
Definition: AliEmcalJetShapeProperties.h:104
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:121
AliEmcalJet::Py
Double_t Py() const
Definition: AliEmcalJet.h:107
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:117
AliAnalysisTaskEmcalQGTagging::GetJetMass
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:621
AliEmcalPythiaInfo::GetPythiaEventWeight
Float_t GetPythiaEventWeight() const
Definition: AliEmcalPythiaInfo.h:56
AliTrackContainer
Container with name, TClonesArray and cuts for particles.
Definition: AliTrackContainer.h:54
AliEmcalJet::GetLabel
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
AliEmcalParticle.h
AliEmcalPythiaInfo::GetPartonPhi7
Float_t GetPartonPhi7() const
Definition: AliEmcalPythiaInfo.h:51
AliAnalysisTaskEmcalQGTagging::GetJetAngularity
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:657
AliEmcalPythiaInfo::GetPartonFlag7
Int_t GetPartonFlag7() const
Definition: AliEmcalPythiaInfo.h:48
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
AliEmcalPythiaInfo::GetPartonPhi6
Float_t GetPartonPhi6() const
Definition: AliEmcalPythiaInfo.h:46
AliEmcalJet::TrackAt
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
AliAnalysisTaskEmcalQGTagging::GetJetCoronna
Float_t GetJetCoronna(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:739
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:63
AliEmcalJet::Px
Double_t Px() const
Definition: AliEmcalJet.h:106
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1623
AliJetContainer::GetRhoMassName
const TString & GetRhoMassName() const
Definition: AliJetContainer.h:168
AliEmcalPythiaInfo::GetPartonFlag6
Int_t GetPartonFlag6() const
Definition: AliEmcalPythiaInfo.h:43
AliJetContainer::GetParticleContainer
AliParticleContainer * GetParticleContainer() const
Definition: AliJetContainer.h:180
AliEmcalJetShapeProperties::GetFirstOrderSubtractedConstituent
Double_t GetFirstOrderSubtractedConstituent() const
Definition: AliEmcalJetShapeProperties.h:103
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskEmcalQGTagging::GetJetCircularity
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:869
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliAnalysisTaskEmcalQGTagging::PTD
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:752
Float_t
float Float_t
Definition: External.C:68
AliParticleContainer.h
AliJetContainer::GetNJets
Int_t GetNJets() const
Definition: AliJetContainer.h:153
AliEmcalJetShapeProperties::GetSecondOrderSubtractedAngularity
Double_t GetSecondOrderSubtractedAngularity() const
Definition: AliEmcalJetShapeProperties.h:63
AliJetContainer.h
AliAnalysisTaskEmcalQGTagging::GetSigma2
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1004
AliParticleContainer::GetAcceptParticle
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Definition: AliParticleContainer.cxx:113
AliAnalysisTaskEmcal::fCent
Double_t fCent
!event centrality
Definition: AliAnalysisTaskEmcal.h:886
AliAnalysisTaskEmcalQGTagging::RecursiveParents
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1108
AliRhoParameter.h
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:299
AliAnalysisTaskEmcalQGTagging::Angularity
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:631
Data
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)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskEmcalQGTagging::LeSub
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:881
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliAnalysisTaskEmcalJet::GetRhoVal
Double_t GetRhoVal(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:626
AliEmcalJetShapeProperties::GetFirstOrderSubtractedCircularity
Double_t GetFirstOrderSubtractedCircularity() const
Definition: AliEmcalJetShapeProperties.h:82
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:909
AliAnalysisTaskEmcalQGTagging.h
AliAnalysisTaskEmcalQGTagging::Sigma2
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:943
AliAnalysisTaskEmcalQGTagging::RelativePhi
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1094
nVar
const Int_t nVar
Definition: AliAnalysisTaskSubJetFraction.h:20
AliEmcalJetShapeProperties::GetSecondOrderSubtractedCircularity
Double_t GetSecondOrderSubtractedCircularity() const
Definition: AliEmcalJetShapeProperties.h:83
AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:930
AliAnalysisTaskEmcal::GetTrackContainer
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcal.h:223
AliEmcalPythiaInfo
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
Definition: AliEmcalPythiaInfo.h:31
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:256
AliTrackContainer::GetAcceptTrack
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Definition: AliTrackContainer.cxx:252
AliAnalysisTaskEmcalJet
Base task in the EMCAL jet framework.
Definition: AliAnalysisTaskEmcalJet.h:30
AliEmcalJet
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
AliAnalysisTaskEmcal::GetPythiaInfo
const AliEmcalPythiaInfo * GetPythiaInfo() const
Definition: AliAnalysisTaskEmcal.h:333
swap
void swap(AliEmcalContainerIndexMap< X, Y > &first, AliEmcalContainerIndexMap< X, Y > &second)
Definition: AliEmcalContainerIndexMap.h:171
AliAnalysisTaskEmcalQGTagging::GetJetpTD
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:775
AliEmcalJet::Pz
Double_t Pz() const
Definition: AliEmcalJet.h:108
AliJetContainer::GetFractionSharedPt
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Definition: AliJetContainer.cxx:806
AliEmcalPythiaInfo.h
Declaration of class AliEmcalPythiaInfo.
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:336
AliEmcalJetShapeProperties::GetFirstOrderSubtractedSigma2
Double_t GetFirstOrderSubtractedSigma2() const
Definition: AliEmcalJetShapeProperties.h:92
AliAnalysisTaskEmcalQGTagging::Circularity
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:786
Bool_t
bool Bool_t
Definition: External.C:53
AliEmcalJet::GetShapeProperties
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
AliAnalysisTaskEmcalQGTagging::GetJetLeSub
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:918
AliEmcalPythiaInfo::GetPartonPt7
Float_t GetPartonPt7() const
Definition: AliEmcalPythiaInfo.h:49
AliEmcalJet::M
Double_t M() const
Definition: AliEmcalJet.h:120
AliEmcalJet.h
AliJetContainer
Container for jet within the EMCAL jet framework.
Definition: AliJetContainer.h:37
AliAnalysisTaskEmcalQGTagging::Run
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Definition: AliAnalysisTaskEmcalQGTagging.cxx:257
TH1
Definition: External.C:196
AliEmcalJet::GetPtSortedTrackConstituentIndexes
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Definition: AliEmcalJet.cxx:424
AliEmcalPythiaInfo::GetPartonPt6
Float_t GetPartonPt6() const
Definition: AliEmcalPythiaInfo.h:44
AliJetContainer::GetJet
AliEmcalJet * GetJet(Int_t i) const
Definition: AliJetContainer.cxx:270