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