AliPhysics  fffcdf3 (fffcdf3)
 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] = "mJet";
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] = "mJetMatch";
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] = GetJetMass(jet1,0);
531  fShapesVar[4] = GetJetAngularity(jet1,0);
532  //fShapesVar[5] = GetJetCircularity(jet1,0);
533  fShapesVar[5] = GetJetLeSub(jet1,0);
534  //fShapesVar[6] = GetJetCoronna(jet1,0);
535  RecursiveParents(jet1,jetCont,0);
536  RecursiveParents(jet1,jetCont,1);
537  RecursiveParents(jet1,jetCont,2);
538 
539  Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0., coronnaMatch=0;
540  Int_t kMatched = 0;
541 
542  if (fJetShapeType==kPythiaDef) {
543  kMatched =1;
544  if(fJetShapeSub==kConstSub) kMatched = 3;
545 
546  ptMatch=jet3->Pt();
547  ptDMatch=GetJetpTD(jet3, kMatched);
548  massMatch=GetJetMass(jet3,kMatched);
549  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
550  angulMatch=GetJetAngularity(jet3, kMatched);
551  //circMatch=GetJetCircularity(jet3, kMatched);
552  lesubMatch=GetJetLeSub(jet3, kMatched);
553  //coronnaMatch=GetJetCoronna(jet3,kMatched);
554  //sigma2Match = GetSigma2(jet2, kMatched);
555  }
556 
557  if (fJetShapeType==kDetEmbPartPythia) {
558  if(fJetShapeSub==kConstSub) kMatched = 3;
559  if(fJetShapeSub==kDerivSub) kMatched = 2;
560  ptMatch=jet3->Pt();
561  ptDMatch=GetJetpTD(jet3, kMatched);
562  massMatch=GetJetMass(jet3,kMatched);
563  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
564  angulMatch=GetJetAngularity(jet3, kMatched);
565  // circMatch=GetJetCircularity(jet3, kMatched);
566  lesubMatch=GetJetLeSub(jet3, kMatched);
567  //coronnaMatch = GetJetCoronna(jet3, kMatched);
568 
569  }
570 
571 
572 
573  if (fJetShapeType == kMCTrue || fJetShapeType == kData || fJetShapeType == kGenOnTheFly) {
574  kMatched = 0;
575  ptMatch=0.;
576  ptDMatch=0.;
577  massMatch=0.;
578  //constMatch=0.;
579  angulMatch=0.;
580  // circMatch=0.;
581  lesubMatch=0.;
582  //coronnaMatch =0.;
583 
584  }
585 
586 
587 
588  fShapesVar[6] = ptMatch;
589  fShapesVar[7] = ptDMatch;
590  fShapesVar[8] = massMatch;
591  fShapesVar[9] = angulMatch;
592  //fShapesVar[12] = circMatch;
593  fShapesVar[10] = lesubMatch;
594  // fShapesVar[14] = coronnaMatch;
595  fShapesVar[11] = kWeight;
596  //fShapesVar[16] = ntracksEvt;
597  // fShapesVar[16] = rhoVal;
598  //fShapesVar[17] = rhoMassVal;
599  //fShapesVar[16] = jet1->Pt();
600 
601 
602  fTreeObservableTagging->Fill();
603 
604 
605 
606 
607 
608 
609  }
610 
611  }
612 
613  return kTRUE;
614 }
615 
616 //________________________________________________________________________
617 Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0){
618  //calc subtracted jet mass
619  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
620  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtracted();
621  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
622  else
623  return jet->M();
624 }
625 
626 //________________________________________________________________________
627 Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
628 
629  AliJetContainer *jetCont = GetJetContainer(jetContNb);
630  if (!jet->GetNumberOfTracks())
631  return 0;
632  Double_t den=0.;
633  Double_t num = 0.;
634  AliVParticle *vp1 = 0x0;
635  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
636  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
637 
638  if (!vp1){
639  Printf("AliVParticle associated to constituent not found");
640  continue;
641  }
642 
643  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
644  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
645  Double_t dr = TMath::Sqrt(dr2);
646  num=num+vp1->Pt()*dr;
647  den=den+vp1->Pt();
648  }
649  return num/den;
650 }
651 
652 //________________________________________________________________________
653 Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0){
654 
655  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
656  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedAngularity();
657  else return jet->GetShapeProperties()->GetSecondOrderSubtractedAngularity();
658  else
659  return Angularity(jet, jetContNb);
660 
661 }
662 
663 //____________________________________________________________________________
664 
665 Float_t AliAnalysisTaskEmcalQGTagging::Coronna(AliEmcalJet *jet, Int_t jetContNb = 0){
666 
667  AliTrackContainer *PartCont = NULL;
668  AliParticleContainer *PartContMC = NULL;
669 
670 
671  if (fJetShapeSub==kConstSub){
672  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(1);
673  else PartCont = GetTrackContainer(1);
674  }
675  else{
676  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(0);
677  else PartCont = GetTrackContainer(0);
678  }
679 
680  TClonesArray *TracksArray = NULL;
681  TClonesArray *TracksArrayMC = NULL;
682 
683  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
684  else TracksArray = PartCont->GetArray();
685 
686  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
687  if(!PartContMC || !TracksArrayMC) return -2;
688  }
689  else {
690  if(!PartCont || !TracksArray) return -2;
691  }
692 
693 
694  AliAODTrack *Track = 0x0;
695  Float_t sumpt=0;
696  Int_t NTracks=0;
697  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
698  else NTracks = TracksArray->GetEntriesFast();
699 
700  for(Int_t i=0; i < NTracks; i++){
701  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
702  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
703  if (!Track) continue;
704  if(TMath::Abs(Track->Eta())>0.9) continue;
705  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
706  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
707  Double_t dr = TMath::Sqrt(dr2);
708  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
709  }
710  }
711  else{
712  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
713  if (!Track) continue;
714  if(TMath::Abs(Track->Eta())>0.9) continue;
715  if (Track->Pt()<0.15) continue;
716  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
717  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
718  Double_t dr = TMath::Sqrt(dr2);
719  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
720 
721  }
722  }
723  }
724 
725 
726 
727  return sumpt;
728 
729 
730 
731 
732 }
733 
734 //________________________________________________________________________
735 Float_t AliAnalysisTaskEmcalQGTagging::GetJetCoronna(AliEmcalJet *jet, Int_t jetContNb = 0){
736 
737  if((fJetShapeSub==kDerivSub) && (jetContNb==0)) return -2;
738  else
739  return Coronna(jet, jetContNb);
740 
741 }
742 
743 
744 
745 
746 
747 //________________________________________________________________________
748 Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
749 
750  AliJetContainer *jetCont = GetJetContainer(jetContNb);
751  if (!jet->GetNumberOfTracks())
752  return 0;
753  Double_t den=0.;
754  Double_t num = 0.;
755  AliVParticle *vp1 = 0x0;
756  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
757  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
758 
759  if (!vp1){
760  Printf("AliVParticle associated to constituent not found");
761  continue;
762  }
763 
764  num=num+vp1->Pt()*vp1->Pt();
765  den=den+vp1->Pt();
766  }
767  return TMath::Sqrt(num)/den;
768 }
769 
770 //________________________________________________________________________
771 Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0){
772  //calc subtracted jet mass
773  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
774  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedpTD();
775  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
776  else
777  return PTD(jet, jetContNb);
778 
779 }
780 
781 //_____________________________________________________________________________
782 Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
783 
784  AliJetContainer *jetCont = GetJetContainer(jetContNb);
785  if (!jet->GetNumberOfTracks())
786  return 0;
787  Double_t mxx = 0.;
788  Double_t myy = 0.;
789  Double_t mxy = 0.;
790  int nc = 0;
791  Double_t sump2 = 0.;
792  Double_t pxjet=jet->Px();
793  Double_t pyjet=jet->Py();
794  Double_t pzjet=jet->Pz();
795 
796 
797  //2 general normalized vectors perpendicular to the jet
798  TVector3 ppJ1(pxjet, pyjet, pzjet);
799  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
800  ppJ3.SetMag(1.);
801  TVector3 ppJ2(-pyjet, pxjet, 0);
802  ppJ2.SetMag(1.);
803  AliVParticle *vp1 = 0x0;
804  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
805  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
806 
807  if (!vp1){
808  Printf("AliVParticle associated to constituent not found");
809  continue;
810  }
811 
812  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
813 
814  //local frame
815  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
816  TVector3 pPerp = pp - pLong;
817  //projection onto the two perpendicular vectors defined above
818 
819  Float_t ppjX = pPerp.Dot(ppJ2);
820  Float_t ppjY = pPerp.Dot(ppJ3);
821  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
822  if(ppjT<=0) return 0;
823 
824  mxx += (ppjX * ppjX / ppjT);
825  myy += (ppjY * ppjY / ppjT);
826  mxy += (ppjX * ppjY / ppjT);
827  nc++;
828  sump2 += ppjT;}
829 
830  if(nc<2) return 0;
831  if(sump2==0) return 0;
832  // Sphericity Matrix
833  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
834  TMatrixDSym m0(2,ele);
835 
836  // Find eigenvectors
837  TMatrixDSymEigen m(m0);
838  TVectorD eval(2);
839  TMatrixD evecm = m.GetEigenVectors();
840  eval = m.GetEigenValues();
841  // Largest eigenvector
842  int jev = 0;
843  // cout<<eval[0]<<" "<<eval[1]<<endl;
844  if (eval[0] < eval[1]) jev = 1;
845  TVectorD evec0(2);
846  // Principle axis
847  evec0 = TMatrixDColumn(evecm, jev);
848  Double_t compx=evec0[0];
849  Double_t compy=evec0[1];
850  TVector2 evec(compx, compy);
851  Double_t circ=0;
852  if(jev==1) circ=2*eval[0];
853  if(jev==0) circ=2*eval[1];
854 
855  return circ;
856 
857 
858 
859 }
860 
861 
862 
863 
864 //________________________________________________________________________
865 Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ){
866  //calc subtracted jet mass
867 
868  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
869  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedCircularity();
870  else return jet->GetShapeProperties()->GetSecondOrderSubtractedCircularity();
871  else
872  return Circularity(jet, jetContNb);
873 
874 }
875 
876 //________________________________________________________________________
877 Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
878 
879  AliJetContainer *jetCont = GetJetContainer(jetContNb);
880  if (!jet->GetNumberOfTracks())
881  return 0;
882  Double_t den=0.;
883  Double_t num = 0.;
884  AliVParticle *vp1 = 0x0;
885  AliVParticle *vp2 = 0x0;
886  std::vector<int> ordindex;
887  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
888  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
889  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
890 
891  if(ordindex.size()<2) return -1;
892 
893  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
894  if (!vp1){
895  Printf("AliVParticle associated to Leading constituent not found");
896  return -1;
897  }
898 
899  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
900  if (!vp2){
901  Printf("AliVParticle associated to Subleading constituent not found");
902  return -1;
903  }
904 
905 
906  num=vp1->Pt();
907  den=vp2->Pt();
908  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
909 
910 return num-den;
911 }
912 
913 //________________________________________________________________________
914 Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
915  //calc subtracted jet mass
916 
917  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
918  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedLeSub();
919  else return jet->GetShapeProperties()->GetSecondOrderSubtractedLeSub();
920  else
921  return LeSub(jet, jetContNb);
922 
923 }
924 
925 //________________________________________________________________________
926  Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0){
927  //calc subtracted jet mass
928 
929  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
930  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedConstituent();
931  else return jet->GetShapeProperties()->GetSecondOrderSubtractedConstituent();
932  else
933  return jet->GetNumberOfTracks();
934 
935  }
936 
937 
938 //______________________________________________________________________________
939 Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
940 
941  AliJetContainer *jetCont = GetJetContainer(jetContNb);
942  if (!jet->GetNumberOfTracks())
943  return 0;
944  Double_t mxx = 0.;
945  Double_t myy = 0.;
946  Double_t mxy = 0.;
947  int nc = 0;
948  Double_t sump2 = 0.;
949 
950  AliVParticle *vp1 = 0x0;
951  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
952  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
953 
954  if (!vp1){
955  Printf("AliVParticle associated to constituent not found");
956  continue;
957  }
958 
959  Double_t ppt=vp1->Pt();
960  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
961 
962  Double_t deta = vp1->Eta()-jet->Eta();
963  mxx += ppt*ppt*deta*deta;
964  myy += ppt*ppt*dphi*dphi;
965  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
966  nc++;
967  sump2 += ppt*ppt;
968 
969  }
970  if(nc<2) return 0;
971  if(sump2==0) return 0;
972  // Sphericity Matrix
973  Double_t ele[4] = {mxx , mxy , mxy , myy };
974  TMatrixDSym m0(2,ele);
975 
976  // Find eigenvectors
977  TMatrixDSymEigen m(m0);
978  TVectorD eval(2);
979  TMatrixD evecm = m.GetEigenVectors();
980  eval = m.GetEigenValues();
981  // Largest eigenvector
982  int jev = 0;
983  // cout<<eval[0]<<" "<<eval[1]<<endl;
984  if (eval[0] < eval[1]) jev = 1;
985  TVectorD evec0(2);
986  // Principle axis
987  evec0 = TMatrixDColumn(evecm, jev);
988  Double_t compx=evec0[0];
989  Double_t compy=evec0[1];
990  TVector2 evec(compx, compy);
991  Double_t sig=0;
992  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
993  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
994 
995  return sig;
996 
997 }
998 
999 //________________________________________________________________________
1000 Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0){
1001  //calc subtracted jet mass
1002 
1003  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1004  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtractedSigma2();
1005  else return jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2();
1006  else
1007  return Sigma2(jet, jetContNb);
1008 
1009 }
1010 
1011 //________________________________________________________________________
1012 Int_t AliAnalysisTaskEmcalQGTagging::SelectTrigger(Float_t minpT, Float_t maxpT){
1013 
1014  AliTrackContainer *PartCont = NULL;
1015  AliParticleContainer *PartContMC = NULL;
1016 
1017 
1018  if (fJetShapeSub==kConstSub){
1019  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(1);
1020  else PartCont = GetTrackContainer(1);
1021  }
1022  else{
1023  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) PartContMC = GetParticleContainer(0);
1024  else PartCont = GetTrackContainer(0);
1025  }
1026 
1027  TClonesArray *TracksArray = NULL;
1028  TClonesArray *TracksArrayMC = NULL;
1029 
1030  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1031  else TracksArray = PartCont->GetArray();
1032 
1033  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
1034  if(!PartContMC || !TracksArrayMC) return -99999;
1035  }
1036  else {
1037  if(!PartCont || !TracksArray) return -99999;
1038  }
1039 
1040 
1041  AliAODTrack *Track = 0x0;
1042 
1043 
1044 
1045  TList *trackList = new TList();
1046  Int_t triggers[100];
1047  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1048  Int_t iTT = 0;
1049  Int_t NTracks=0;
1050  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1051  else NTracks = TracksArray->GetEntriesFast();
1052 
1053  for(Int_t i=0; i < NTracks; i++){
1054  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly){
1055  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1056  if (!Track) continue;
1057  if(TMath::Abs(Track->Eta())>0.9) continue;
1058  if (Track->Pt()<0.15) continue;
1059  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1060  triggers[iTT] = i;
1061  iTT++;
1062  }
1063  }
1064  }
1065  else{
1066  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1067  if (!Track) continue;
1068  if(TMath::Abs(Track->Eta())>0.9) continue;
1069  if (Track->Pt()<0.15) continue;
1070  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1071  triggers[iTT] = i;
1072  iTT++;
1073  }
1074  }
1075  }
1076  }
1077 
1078 
1079  if (iTT == 0) return -99999;
1080  Int_t nbRn = 0, index = 0 ;
1081  TRandom3* random = new TRandom3(0);
1082  nbRn = random->Integer(iTT);
1083  index = triggers[nbRn];
1084  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1085  return index;
1086 
1087 }
1088 
1089 //__________________________________________________________________________________
1090 Double_t AliAnalysisTaskEmcalQGTagging::RelativePhi(Double_t mphi,Double_t vphi){
1091 
1092  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1093  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1094  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1095  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1096  double dphi = mphi-vphi;
1097  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1098  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1099  return dphi;//dphi in [-Pi, Pi]
1100 }
1101 
1102 
1103 //_________________________________________________________________________
1104 void AliAnalysisTaskEmcalQGTagging::RecursiveParents(AliEmcalJet *fJet,AliJetContainer *fJetCont, Int_t ReclusterAlgo){
1105 
1106  std::vector<fastjet::PseudoJet> fInputVectors;
1107  fInputVectors.clear();
1108  fastjet::PseudoJet PseudoTracks;
1109  double xflagalgo=0;
1110  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1111 
1112  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1113  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1114  if (!fTrk) continue;
1115  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1116  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1117  fInputVectors.push_back(PseudoTracks);
1118 
1119  }
1120  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1121 
1122  if(ReclusterAlgo==0){ xflagalgo=0.5;
1123  jetalgo=fastjet::kt_algorithm ;}
1124 
1125  if(ReclusterAlgo==1){ xflagalgo=1.5;
1126  jetalgo=fastjet::cambridge_algorithm;}
1127  if(ReclusterAlgo==2){ xflagalgo=2.5;
1128  jetalgo=fastjet::antikt_algorithm;}
1129 
1130  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1131 
1132  try {
1133  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1134  std::vector<fastjet::PseudoJet> fOutputJets;
1135  fOutputJets.clear();
1136  fOutputJets=fClustSeqSA.inclusive_jets(0);
1137 
1138  fastjet::PseudoJet jj;
1139  fastjet::PseudoJet j1;
1140  fastjet::PseudoJet j2;
1141  jj=fOutputJets[0];
1142  double ndepth=0;
1143  while(jj.has_parents(j1,j2)){
1144  ndepth=ndepth+1;
1145  if(j1.perp() < j2.perp()) swap(j1,j2);
1146  double delta_R=j1.delta_R(j2);
1147  double z=j2.perp()/(j1.perp()+j2.perp());
1148  double y =log(1.0/delta_R);
1149  double lnpt_rel=log(z*delta_R);
1150  Double_t LundEntries[5] = {y,lnpt_rel,fOutputJets[0].perp(),xflagalgo,ndepth};
1151  fHLundIterative->Fill(LundEntries);
1152  jj=j1;}
1153 
1154 
1155 
1156 
1157  } catch (fastjet::Error) {
1158  AliError(" [w] FJ Exception caught.");
1159  //return -1;
1160  }
1161 
1162 
1163 
1164 
1165  return;
1166 
1167 
1168 }
1169 
1170 
1171 
1172 
1173 
1174 
1175 
1176 //________________________________________________________________________
1177 Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
1178  //
1179  // retrieve event objects
1180  //
1181  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
1182  return kFALSE;
1183 
1184  return kTRUE;
1185 }
1186 
1187 //_______________________________________________________________________
1188 void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *)
1189 {
1190  // Called once at the end of the analysis.
1191 
1192  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
1193  // if (!fTreeObservableTagging){
1194  // Printf("ERROR: fTreeObservableTagging not available");
1195  // return;
1196  // }
1197 
1198 }
1199 
AliAnalysisTaskEmcalQGTagging::SelectTrigger
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1012
AliAnalysisTaskEmcalQGTagging::Coronna
Float_t Coronna(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:665
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:1177
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:617
AliEmcalPythiaInfo::GetPythiaEventWeight
Float_t GetPythiaEventWeight() const
Definition: AliEmcalPythiaInfo.h:56
AliTrackContainer
Container with name, TClonesArray and cuts for particles.
Definition: AliTrackContainer.h:29
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:653
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:735
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:865
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliAnalysisTaskEmcalQGTagging::PTD
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:748
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:1000
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:1104
AliRhoParameter.h
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:299
AliAnalysisTaskEmcalQGTagging::Angularity
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:627
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:877
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:939
AliAnalysisTaskEmcalQGTagging::RelativePhi
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Definition: AliAnalysisTaskEmcalQGTagging.cxx:1090
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:926
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:210
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:771
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:782
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:914
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