AliPhysics  9c7580a (9c7580a)
AliAnalysisTaskEmcalJetShapesMC.cxx
Go to the documentation of this file.
1 //
2 // Jet QG tagging analysis task.
3 //
4 // Author: D. Caffarri, L. Cunqueiro
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TTree.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TProfile.h>
15 #include <TChain.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TKey.h>
19 #include <AliAnalysisDataSlot.h>
20 #include <AliAnalysisDataContainer.h>
21 #include "TMatrixD.h"
22 #include "TMatrixDSym.h"
23 #include "TMatrixDSymEigen.h"
24 #include "TVector3.h"
25 #include "TVector2.h"
26 #include "AliVCluster.h"
27 #include "AliVTrack.h"
28 #include "AliEmcalJet.h"
29 #include "AliRhoParameter.h"
30 #include "AliLog.h"
31 #include "AliEmcalParticle.h"
32 #include "AliMCEvent.h"
33 #include "AliGenPythiaEventHeader.h"
34 #include "AliAODMCHeader.h"
35 #include "AliMCEvent.h"
36 #include "AliAnalysisManager.h"
37 #include "AliJetContainer.h"
38 #include "AliParticleContainer.h"
39 #include "AliEmcalPythiaInfo.h"
40 #include "TRandom3.h"
41 #include "TF1.h"
42 #include "AliEmcalJetFinder.h"
43 #include "AliAODEvent.h"
45 
46 using std::cout;
47 using std::endl;
48 
50 
51 //________________________________________________________________________
54  fContainer(0),
55  fMinFractionShared(0),
56  fJetShapeType(kGenShapes),
57  fJetShapeSub(kNoSub),
58  fJetSelection(kInclusive),
59  fPtThreshold(-9999.),
60  fRMatching(0.2),
61  fJetRadius(0.4),
62  fSubjetRadius(0.2),
63  fSelectedShapes(0),
64  fSwitchKtNSub(0),
65  fSwitchMinNSub(0),
66  fSwitchAktNSub(0),
67  fSwitchSDKtNSub(0),
68  fSwitchSDMinNSub(0),
69  fAdditionalTracks(0),
70  fHardCutoff(0),
71  fminpTTrig(20.),
72  fmaxpTTrig(50.),
73  fangWindowRecoil(0.6),
74  fSemigoodCorrect(0),
75  fHolePos(0),
76  fHoleWidth(0),
77  fRandom(0),
78  fqhat(1),
79  fxlength(2),
80  fCentSelectOn(kTRUE),
81  fCentMin(0),
82  fCentMax(10),
83  fOneConstSelectOn(kFALSE),
84  fDerivSubtrOrder(0),
85  fPhiJetCorr6(0x0),
86  fPhiJetCorr7(0x0),
87  fEtaJetCorr6(0x0),
88  fEtaJetCorr7(0x0),
89  fPtJetCorr(0x0),
90  fPtJet(0x0),
91  fhpTjetpT(0x0),
92  fhPt(0x0),
93  fhPhi(0x0),
94  fHLundIterative(0x0),
95  fHLundIterativeInject(0x0),
96  fNbOfConstvspT(0x0),
97  fTreeObservableTagging(0x0),
98  fTf1Omega(0x0),
99  fTf1Kt(0x0),
100  fScaleELoss(kFALSE),
101  xfraction(1),
102  fAddMedScat(kFALSE),
103  fAddMedScatPtFrac(1),
104  fAddMedScatN(100)
105 
106 
107 {
108  for(Int_t i=0;i<11;i++){
109  fShapesVar[i]=0;}
110  SetMakeGeneralHistograms(kTRUE);
111  DefineOutput(1, TList::Class());
112  DefineOutput(2, TTree::Class());
113 }
114 
115 //________________________________________________________________________
117  AliAnalysisTaskEmcalJet(name, kTRUE),
118  fContainer(0),
119  fMinFractionShared(0),
120  fJetShapeType(kGenShapes),
121  fJetShapeSub(kNoSub),
122  fJetSelection(kInclusive),
123  fPtThreshold(-9999.),
124  fRMatching(0.2),
125  fSelectedShapes(0),
126  fSwitchKtNSub(0),
127  fSwitchMinNSub(0),
128  fSwitchAktNSub(0),
129  fSwitchSDKtNSub(0),
130  fSwitchSDMinNSub(0),
131  fAdditionalTracks(0),
132  fHardCutoff(0),
133  fminpTTrig(20.),
134  fmaxpTTrig(50.),
135  fangWindowRecoil(0.6),
136  fSemigoodCorrect(0),
137  fHolePos(0),
138  fHoleWidth(0),
139  fRandom(0),
140  fqhat(1),
141  fxlength(2),
142  fCentSelectOn(kTRUE),
143  fCentMin(0),
144  fCentMax(10),
145  fOneConstSelectOn(kFALSE),
146  fDerivSubtrOrder(0),
147  fPhiJetCorr6(0x0),
148  fPhiJetCorr7(0x0),
149  fEtaJetCorr6(0x0),
150  fEtaJetCorr7(0x0),
151  fPtJetCorr(0x0),
152  fPtJet(0x0),
153  fhpTjetpT(0x0),
154  fhPt(0x0),
155  fhPhi(0x0),
156  fHLundIterative(0x0),
157  fHLundIterativeInject(0x0),
158  fNbOfConstvspT(0x0),
159  fTreeObservableTagging(0x0),
160  fTf1Omega(0x0),
161  fTf1Kt(0x0),
162  fScaleELoss(kFALSE),
163  xfraction(1),
164  fAddMedScat(kFALSE),
165  fAddMedScatPtFrac(1),
166  fAddMedScatN(100)
167 {
168  // Standard constructor.
169 
170 
171  for(Int_t i=0;i<11;i++){
172  fShapesVar[i]=0;}
173 
175 
176  DefineOutput(1, TList::Class());
177  DefineOutput(2, TTree::Class());
178 
179 
180 }
181 
182 //________________________________________________________________________
184 {
186  delete fTreeObservableTagging;
188  }
189 
190  if(fRandom) delete fRandom;
191  if(fTf1Omega) delete fTf1Omega;
192  if(fTf1Kt) delete fTf1Kt;
193 
194 }
195 
196 //________________________________________________________________________
198 {
199  // Create user output.
201 
202  Bool_t oldStatus = TH1::AddDirectoryStatus();
203  TH1::AddDirectory(oldStatus);
204 
205  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
206 
207  //TH1::AddDirectory(oldStatus);
208 
209  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
210  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
211 
212  const Int_t nVar = 11;
213 
214  TString *fShapesVarNames = new TString [nVar];
215 
216  fShapesVarNames[0] = "partonCode";
217  fShapesVarNames[1] = "ptJet";
218  fShapesVarNames[2] = "ptDJet";
219  fShapesVarNames[3] = "mJet";
220  fShapesVarNames[4] = "nbOfConst";
221  fShapesVarNames[5] = "angularity";
222  fShapesVarNames[6] = "nitersd";
223  fShapesVarNames[7] = "niterall";
224  fShapesVarNames[8] = "weightPythia";
225  //fShapesVarNames[6] = "Nsubjet1kt";
226  // fShapesVarNames[7] = "Nsubjet2kt";
227  // fShapesVarNames[8] = "Nsubjet1Min";
228  // fShapesVarNames[9] = "Nsubjet2Min";
229  // fShapesVarNames[10] = "DeltaRkt";
230  // fShapesVarNames[11] = "DeltaRMin";
231  fShapesVarNames[9] = "SDSymm";
232  fShapesVarNames[10] = "scaledptJet";
233  // fShapesVarNames[13] = "SDDeltaR";
234  // fShapesVarNames[14] = "SDGroomedFrac";
235  // fShapesVarNames[15] = "SDGroomedN";
236  // fShapesVarNames[16] = "SDMass";
237  // fShapesVarNames[17] = "SDSymmkt";
238  // fShapesVarNames[18] = "SDDeltaRkt";
239  // fShapesVarNames[19] = "SDGroomedFrackt";
240  // fShapesVarNames[20] = "SDGroomedNkt";
241  // fShapesVarNames[21] = "SDMasskt";
242  // fShapesVarNames[22] = "SDSymmAkt";
243  // fShapesVarNames[23] = "SDDeltaRAkt";
244  // fShapesVarNames[24] = "SDGroomedFracAkt";
245  // fShapesVarNames[25] = "SDGroomedNAkt";
246  // fShapesVarNames[26] = "SDMassAkt";
247  // fShapesVarNames[27] = "SDSymmktForm";
248  // fShapesVarNames[28] = "SDDeltaRktForm";
249  // fShapesVarNames[29] = "SDGroomedFracktForm";
250  // fShapesVarNames[30] = "SDGroomedNktForm";
251  // fShapesVarNames[31] = "SDMassktForm";
252  // fShapesVarNames[32] = "SDSymmDemo";
253  // fShapesVarNames[33] = "SDDeltaRDemo";
254  // fShapesVarNames[34] = "SDGroomedFracDemo";
255  // fShapesVarNames[35] = "SDGroomedNDemo";
256  // fShapesVarNames[36] = "SDMassDemo";
257  // fShapesVarNames[42] = "SDSymmForm";
258  // fShapesVarNames[43] = "SDDeltaRForm";
259  // fShapesVarNames[44] = "SDGroomedFracForm";
260  // fShapesVarNames[45] = "SDGroomedNForm";
261  // fShapesVarNames[46] = "SDMassForm";
262  // fShapesVarNames[47] = "weightPythia";
263  // fShapesVarNames[42] = "SDSymmNoCut";
264  // fShapesVarNames[43] = "SDDeltaRNoCut";
265  // fShapesVarNames[44] = "SDGroomedFracNoCut";
266  // fShapesVarNames[45] = "SDGroomedNNoCut";
267  // fShapesVarNames[46] = "SDMassNoCut";
268 
269 
270  //fShapesVarNames[7] = "lesub";
271  //fShapesVarNames[8] = "CoreFraction";
272  //fShapesVarNames[9] = "Nsubjet1";
273  //fShapesVarNames[10] = "Nsubjet2";
274  //fShapesVarNames[11] = "DeltaR";
275  //fShapesVarNames[12] = "OpenAngle";
276  //fShapesVarNames[13] = "weightPythia";
277 
278  //fShapesVarNames[14] = "NT70";
279  //fShapesVarNames[15] = "nConstNT70";
280  //fShapesVarNames[16] = "NT80";
281  //fShapesVarNames[17] = "nConstNT80";
282  //fShapesVarNames[18] = "NT90";
283  //fShapesVarNames[19] = "nConstNT90";
284  //fShapesVarNames[20] = "NT95";
285  //fShapesVarNames[21] = "nConstNT95";
286 
287  //fShapesVarNames[22] = "SubjetFraction";
288 
289 
290  for(Int_t ivar=0; ivar < nVar; ivar++){
291  cout<<"looping over variables"<<endl;
292  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
293 
294 
295 
296  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
297  fOutput->Add(fPhiJetCorr6);
298  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
299  fOutput->Add(fEtaJetCorr6);
300 
301  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
302  fOutput->Add(fPhiJetCorr7);
303  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
304  fOutput->Add(fEtaJetCorr7);
305 
306  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
307  fOutput->Add(fPtJetCorr);
308  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
309  fOutput->Add(fPtJet);
310 
311  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
312  fOutput->Add(fhpTjetpT);
313  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
314  fOutput->Add(fhPt);
315  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
316  fOutput->Add(fhPhi);
317 
318 
319 
320  //log(1/theta),log(z*theta),jetpT,algo//
321  const Int_t dimSpec = 6;
322  const Int_t nBinsSpec[6] = {50,50,10,3,22,10};
323  const Double_t lowBinSpec[6] = {0.0,-10, 0,0,0,0};
324  const Double_t hiBinSpec[6] = {5.0, 0,200,3,22,10};
325  fHLundIterative = new THnSparseF("fHLundIterative",
326  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,partonFlavor,depth]",
327  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
328  fOutput->Add(fHLundIterative);
329 
330  if(fAdditionalTracks>0){
331  //log(1/theta),log(z*theta),jetpT of added tracks
332  const Int_t dimSpecb = 4;
333  const Int_t nBinsSpecb[4] = {50,50,20,20};
334  const Double_t lowBinSpecb[4] = {0.0,-10, 0,0};
335  const Double_t hiBinSpecb[4] = {5.0, 0,200,200};
336  fHLundIterativeInject = new THnSparseF("fHLundIterativeInject",
337  "LundIterativePlotInject [log(1/theta),log(z*theta),pTjet,algo]",
338  dimSpecb,nBinsSpecb,lowBinSpecb,hiBinSpecb);
340 
341 
342  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
343  fOutput->Add(fNbOfConstvspT);
344 
345  //fOutput->Add(fTreeObservableTagging);
346 
347  fRandom = new TRandom3(0);
348  PostData(1, fOutput); // Post data for ALL output slots > 0 here
349  PostData(2, fTreeObservableTagging);
350 
351  delete [] fShapesVarNames;
352 
353  }
354 
355 //________________________________________________________________________
357 {
358  // Run analysis code here, if needed. It will be executed before FillHistograms().
359  // if (gRandom) delete gRandom;
360  // gRandom = new TRandom3(0);
361  return kTRUE;
362 }
363 
364 //________________________________________________________________________
366 {
367  // Fill histograms.
368  //cout<<"IntoFillHistograms"<<endl;
369  AliEmcalJet* jet1 = NULL;
370  AliJetContainer *jetCont = GetJetContainer(0);
371 
372  Float_t kWeight=1;
373  if (fCentSelectOn)
374  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
375 
376  AliAODTrack *triggerHadron = 0x0;
377 
378  if (fJetSelection == kRecoil) {
379  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
380  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
381 
382 
383  if (triggerHadronLabel==-99999) {
384  //Printf ("Trigger Hadron not found, return");
385  return 0;}
386 
387 
389  TClonesArray *trackArrayAn = partContAn->GetArray();
390  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
391 
392  if (!triggerHadron) {
393  //Printf("No Trigger hadron with the found label!!");
394  return 0;
395  }
396 
397  if(fSemigoodCorrect){
398  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
399  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
400  return 0;}
401  }
402 
403  fhPt->Fill(triggerHadron->Pt());
404 
405  }
406 
407  if(jetCont) {
408  jetCont->ResetCurrentID();
409  while((jet1 = jetCont->GetNextAcceptJet())) {
410  //Printf("jet1=%p", jet1);
411  if (!jet1) continue;
412 
413  fPtJet->Fill(jet1->Pt());
414 
415 
416 
417 
419  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
420  if(TMath::Abs(disthole)<fHoleWidth){
421  continue;
422  }
423  }
424 
425  Float_t dphiRecoil = 0.;
426  if (fJetSelection == kRecoil){
427  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
428  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
429  // Printf("Recoil jets back to back not found! continuing");
430  continue;
431  }
432 
433  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
434  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
435  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
436 
437  }
438 
439 
440  fShapesVar[0] = 0.;
441 
442  if (fJetShapeType == kGenShapes){
443  const AliEmcalPythiaInfo *partonsInfo = 0x0;
444  partonsInfo = GetPythiaInfo();
445  //Printf("partonsInfo=%p", partonsInfo);
446  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
447  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
448  kWeight=partonsInfo->GetPythiaEventWeight();
449  //Printf("kWeight=%f", kWeight);
450  fShapesVar[8] = kWeight;
451 
452  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
453  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
454  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
455  if(dRp1 < fRMatching) {
456  fShapesVar[0] = partonsInfo->GetPartonFlag6();
457  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
458  }
459  else {
460  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
461  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
462  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
463  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
464  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
465  if(dRp1 < fRMatching) {
466  fShapesVar[0] = partonsInfo->GetPartonFlag7();
467  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
468  }
469  else fShapesVar[0]=0;
470  }
471  }
472 
473  Double_t ptSubtracted = 0;
474  ptSubtracted= jet1->Pt();
475  //Printf("ptSubtracted=%f", ptSubtracted);
476 
477 
478  if (ptSubtracted < fPtThreshold) continue;
479 
480  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() < 2)) continue;
481 
482  //AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
483  //Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
484 
485 
486 
487 
488 
489 
490 
491  fShapesVar[1] = ptSubtracted;
492  fShapesVar[2] = GetJetpTD(jet1,0);
493  fShapesVar[3] = GetJetMass(jet1,0);
494  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
495  fShapesVar[5] = GetJetAngularity(jet1,0);
496 
497 
498  RecursiveParents(jet1,jetCont,0,fShapesVar[0]);
499  RecursiveParents(jet1,jetCont,1,fShapesVar[0]);
500  RecursiveParents(jet1,jetCont,2,fShapesVar[0]);
501 
502 
503  fTreeObservableTagging->Fill();
504 
505  }
506 
507  }
508 
509  return kTRUE;
510 }
511 
512 //________________________________________________________________________
514  //calc subtracted jet mass
515  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
517  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
518  else
519  return jet->M();
520 }
521 
522 //________________________________________________________________________
524 
525  AliJetContainer *jetCont = GetJetContainer(jetContNb);
526  if (!jet->GetNumberOfTracks())
527  return 0;
528  Double_t den=0.;
529  Double_t num = 0.;
530  AliVParticle *vp1 = 0x0;
531  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
532  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
533 
534  if (!vp1){
535  Printf("AliVParticle associated to constituent not found");
536  continue;
537  }
538 
539  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
540  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
541  Double_t dr = TMath::Sqrt(dr2);
542  num=num+vp1->Pt()*dr;
543  den=den+vp1->Pt();
544  }
545  return num/den;
546 }
547 
548 //________________________________________________________________________
550 
551  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
554  else
555  return Angularity(jet, jetContNb);
556 
557 }
558 
559 
560 //________________________________________________________________________
562 
563  AliJetContainer *jetCont = GetJetContainer(jetContNb);
564  if (!jet->GetNumberOfTracks())
565  return 0;
566  Double_t den=0.;
567  Double_t num = 0.;
568  AliVParticle *vp1 = 0x0;
569  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
570  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
571 
572  if (!vp1){
573  Printf("AliVParticle associated to constituent not found");
574  continue;
575  }
576 
577  num=num+vp1->Pt()*vp1->Pt();
578  den=den+vp1->Pt();
579  }
580  return TMath::Sqrt(num)/den;
581 }
582 
583 //________________________________________________________________________
585  //calc subtracted jet mass
586  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
588  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
589  else
590  return PTD(jet, jetContNb);
591 
592 }
593 
594 //_____________________________________________________________________________
596 
597  AliJetContainer *jetCont = GetJetContainer(jetContNb);
598  if (!jet->GetNumberOfTracks())
599  return 0;
600  Double_t mxx = 0.;
601  Double_t myy = 0.;
602  Double_t mxy = 0.;
603  int nc = 0;
604  Double_t sump2 = 0.;
605  Double_t pxjet=jet->Px();
606  Double_t pyjet=jet->Py();
607  Double_t pzjet=jet->Pz();
608 
609 
610  //2 general normalized vectors perpendicular to the jet
611  TVector3 ppJ1(pxjet, pyjet, pzjet);
612  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
613  ppJ3.SetMag(1.);
614  TVector3 ppJ2(-pyjet, pxjet, 0);
615  ppJ2.SetMag(1.);
616  AliVParticle *vp1 = 0x0;
617  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
618  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
619 
620  if (!vp1){
621  Printf("AliVParticle associated to constituent not found");
622  continue;
623  }
624 
625  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
626 
627  //local frame
628  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
629  TVector3 pPerp = pp - pLong;
630  //projection onto the two perpendicular vectors defined above
631 
632  Float_t ppjX = pPerp.Dot(ppJ2);
633  Float_t ppjY = pPerp.Dot(ppJ3);
634  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
635  if(ppjT<=0) return 0;
636 
637  mxx += (ppjX * ppjX / ppjT);
638  myy += (ppjY * ppjY / ppjT);
639  mxy += (ppjX * ppjY / ppjT);
640  nc++;
641  sump2 += ppjT;}
642 
643  if(nc<2) return 0;
644  if(sump2==0) return 0;
645  // Sphericity Matrix
646  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
647  TMatrixDSym m0(2,ele);
648 
649  // Find eigenvectors
650  TMatrixDSymEigen m(m0);
651  TVectorD eval(2);
652  TMatrixD evecm = m.GetEigenVectors();
653  eval = m.GetEigenValues();
654  // Largest eigenvector
655  int jev = 0;
656  // cout<<eval[0]<<" "<<eval[1]<<endl;
657  if (eval[0] < eval[1]) jev = 1;
658  TVectorD evec0(2);
659  // Principle axis
660  evec0 = TMatrixDColumn(evecm, jev);
661  Double_t compx=evec0[0];
662  Double_t compy=evec0[1];
663  TVector2 evec(compx, compy);
664  Double_t circ=0;
665  if(jev==1) circ=2*eval[0];
666  if(jev==0) circ=2*eval[1];
667 
668  return circ;
669 
670 
671 
672 }
673 
674 
675 
676 
677 //________________________________________________________________________
679  //calc subtracted jet mass
680 
681  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
684  else
685  return Circularity(jet, jetContNb);
686 
687 }
688 
689 //________________________________________________________________________
691 
692  AliJetContainer *jetCont = GetJetContainer(jetContNb);
693  if (!jet->GetNumberOfTracks())
694  return 0;
695  Double_t den=0.;
696  Double_t num = 0.;
697  AliVParticle *vp1 = 0x0;
698  AliVParticle *vp2 = 0x0;
699  std::vector<int> ordindex;
700  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
701  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
702  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
703 
704  if(ordindex.size()<2) return -1;
705 
706  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
707  if (!vp1){
708  Printf("AliVParticle associated to Leading constituent not found");
709  return -1;
710  }
711 
712  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
713  if (!vp2){
714  Printf("AliVParticle associated to Subleading constituent not found");
715  return -1;
716  }
717 
718 
719  num=vp1->Pt();
720  den=vp2->Pt();
721  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
722 
723 return num-den;
724 }
725 
726 //________________________________________________________________________
728  //calc subtracted jet mass
729 
730  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
733  else
734  return LeSub(jet, jetContNb);
735 
736 }
737 
738 //________________________________________________________________________
740  //calc subtracted jet mass
741 
742  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
745  else
746  return jet->GetNumberOfTracks();
747 
748  }
749 
750 
751 //________________________________________________________________________
752 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
753 
754  AliJetContainer *JetCont = GetJetContainer(JetContNb);
755  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
756  Reclusterer->SetRadius(SubJetRadius);
757  Reclusterer->SetJetMinPt(SubJetMinPt);
758  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
759  Reclusterer->SetJetMaxEta(0.9-JetRadius);
760  Reclusterer->SetRecombSheme(0);
761 
762 
763  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
764  Double_t dVtx[3]={0.,0.,0.};
765  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
766  return Reclusterer;
767 }
768 
769 
770 
771 
772 //----------------------------------------------------------------------
774  AliJetContainer *JetCont = GetJetContainer(JetContNb);
775 
776 
777 
778  Double_t SubJetiness_Numerator = 0;
779  Double_t SubJetiness_Denominator = 0;
780  Double_t Index=-2;
781  if (Reclusterer->GetNumberOfJets() < 1) return -2;
782  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
783  if(Index==-999) return -2;
784  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
785  SubJetiness_Denominator=Jet->Pt();
786  return SubJetiness_Numerator/SubJetiness_Denominator;
787 
788 
789 }
790 //__________________________________________________________________________________
792 
793  AliJetContainer *jetCont = GetJetContainer(jetContNb);
794  if (!jet->GetNumberOfTracks())
795  return 0;
796  Double_t den=0.;
797  Double_t num = 0.;
798  AliVParticle *vp1 = 0x0;
799  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
800  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
801 
802  if (!vp1){
803  Printf("AliVParticle associated to constituent not found");
804  continue;
805  }
806 
807  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
808  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
809  Double_t dr = TMath::Sqrt(dr2);
810  if(dr<=fSubjetRadius) num=num+vp1->Pt();
811 
812  }
813  return num/jet->Pt();
814 }
815 
816 
817 
818 
819 //________________________________________________________________________
821  //calc subtracted jet mass
822 
823  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
826  else
827  return CoreFrac(jet, jetContNb);
828 
829 }
830 
831 
832 
833 
834 //----------------------------------------------------------------------
836  AliEmcalJet *SubJet=NULL;
837  Double_t SortingVariable;
838  Int_t ArraySize =N+1;
839  TArrayD *JetSorter = new TArrayD(ArraySize);
840  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
841  for (Int_t i=0; i<ArraySize; i++){
842  JetSorter->SetAt(0,i);
843  }
844  for (Int_t i=0; i<ArraySize; i++){
845  JetIndexSorter->SetAt(0,i);
846  }
847  if(Reclusterer->GetNumberOfJets()<N) return -999;
848  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
849  SubJet=Reclusterer->GetJet(i);
850  if (Type==0) SortingVariable=SubJet->Pt();
851  else if (Type==1) SortingVariable=SubJet->E();
852  else if (Type==2) SortingVariable=SubJet->M();
853  for (Int_t j=0; j<N; j++){
854  if (SortingVariable>JetSorter->GetAt(j)){
855  for (Int_t k=N-1; k>=j; k--){
856  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
857  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
858  }
859  JetSorter->SetAt(SortingVariable,j);
860  JetIndexSorter->SetAt(i,j);
861  break;
862  }
863  }
864  }
865  if (!Index) return JetSorter->GetAt(N-1);
866  else return JetIndexSorter->GetAt(N-1);
867 }
868 
869 
870 
871 //returns -1 if the Nth hardest jet is requested where N>number of available jets
872 //type: 0=Pt 1=E 2=M
873 //Index TRUE=returns index FALSE=returns value of quantatiy in question
874 
875 
876 
877 
878 
879 //______________________________________________________________________________
881 
882  AliJetContainer *jetCont = GetJetContainer(jetContNb);
883  if (!jet->GetNumberOfTracks())
884  return 0;
885  Double_t mxx = 0.;
886  Double_t myy = 0.;
887  Double_t mxy = 0.;
888  int nc = 0;
889  Double_t sump2 = 0.;
890 
891  AliVParticle *vp1 = 0x0;
892  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
893  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
894 
895  if (!vp1){
896  Printf("AliVParticle associated to constituent not found");
897  continue;
898  }
899 
900  Double_t ppt=vp1->Pt();
901  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
902 
903  Double_t deta = vp1->Eta()-jet->Eta();
904  mxx += ppt*ppt*deta*deta;
905  myy += ppt*ppt*dphi*dphi;
906  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
907  nc++;
908  sump2 += ppt*ppt;
909 
910  }
911  if(nc<2) return 0;
912  if(sump2==0) return 0;
913  // Sphericity Matrix
914  Double_t ele[4] = {mxx , mxy , mxy , myy };
915  TMatrixDSym m0(2,ele);
916 
917  // Find eigenvectors
918  TMatrixDSymEigen m(m0);
919  TVectorD eval(2);
920  TMatrixD evecm = m.GetEigenVectors();
921  eval = m.GetEigenValues();
922  // Largest eigenvector
923  int jev = 0;
924  // cout<<eval[0]<<" "<<eval[1]<<endl;
925  if (eval[0] < eval[1]) jev = 1;
926  TVectorD evec0(2);
927  // Principle axis
928  evec0 = TMatrixDColumn(evecm, jev);
929  Double_t compx=evec0[0];
930  Double_t compy=evec0[1];
931  TVector2 evec(compx, compy);
932  Double_t sig=0;
933  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
934  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
935 
936  return sig;
937 
938 }
939 
940 //________________________________________________________________________
942  //calc subtracted jet mass
943 
944  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
947  else
948  return Sigma2(jet, jetContNb);
949 
950 }
951 
952 
953 //_________________________________________________________________________
955 
956  AliJetContainer *jetCont = GetJetContainer(jetContNb);
957  if (!jet->GetNumberOfTracks())
958  return;
959 
960  Double_t ptJet = jet->Pt();
961 
962  AliVParticle *vp1 = 0x0;
963  std::vector<int> ordindex;
964  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
965  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
966  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
967  //if(ordindex.size()<2) return -1;
968 
969  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
970 
971  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
972  if (!vp1){
973  Printf("AliVParticle associated to Leading constituent not found");
974  return;
975  }
976 
977  if (nTFractions[0] <= 0.7*ptJet){
978  nTFractions[0] += vp1->Pt();
979  nTFractions[1] +=1;
980  }
981 
982  if (nTFractions[2] <= 0.8*ptJet){
983  nTFractions[2] += vp1->Pt();
984  nTFractions[3] +=1;
985  }
986 
987  if (nTFractions[4] <= 0.9*ptJet){
988  nTFractions[4] += vp1->Pt();
989  nTFractions[5] +=1;
990  }
991 
992  if (nTFractions[6] <= 0.95*ptJet){
993  nTFractions[6] += vp1->Pt();
994  nTFractions[7] +=1;
995  }
996  }
997 }
998 //_________________________________________________________________________________________________
1000 
1001  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1002 
1003  //Algorithm==0 -> kt_axes;
1004  // Algorithm==1 -> ca_axes;
1005  //Algorithm==2 -> antikt_0p2_axes;
1006  //Algorithm==3 -> wta_kt_axes;
1007  //Algorithm==4 -> wta_ca_axes;
1008  //Algorithm==5 -> onepass_kt_axes;
1009  //Algorithm==6 -> onepass_ca_axes;
1010  //Algorithm==7 -> onepass_antikt_0p2_axes;
1011  //Algorithm==8 -> onepass_wta_kt_axes;
1012  //Algorithm==9 -> onepass_wta_ca_axes;
1013  //Algorithm==10 -> min_axes;
1014 
1015 
1016  //Option==0 returns Nsubjettiness Value
1017  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1018  //Option==2 && N==2 returns Delta R
1019 
1020  if (Jet->GetNumberOfTracks()>=N){
1021  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1022  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1023  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1024  JetFinder->SetRadius(fJetRadius);
1025  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1026  JetFinder->SetRecombSheme(0);
1027  JetFinder->SetJetMinPt(Jet->Pt());
1028  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1029  Double_t dVtx[3]={1,1,1};
1030  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
1031  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,0.2,Beta,Option,0,Beta_SD,ZCut,SoftDropOn);
1032 
1033  }
1034  else return -2;
1035 }
1036 
1037 
1038 
1039 //________________________________________________________________________
1041 
1043  TClonesArray *tracksArray = partCont->GetArray();
1044 
1045  if(!partCont || !tracksArray) return -99999;
1046  AliAODTrack *track = 0x0;
1047  AliEmcalParticle *emcPart = 0x0;
1048 
1049 
1050  TList *trackList = new TList();
1051  Int_t triggers[100];
1052  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1053  Int_t iTT = 0;
1054 
1055  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
1056 
1057 
1058  if (fJetShapeSub == kConstSub){
1059  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
1060  if (!emcPart) continue;
1061  if(TMath::Abs(emcPart->Eta())>0.9) continue;
1062  if (emcPart->Pt()<0.15) continue;
1063 
1064  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
1065  trackList->Add(emcPart);
1066  triggers[iTT] = iTrack;
1067  iTT++;
1068  }
1069  }
1070  else{
1071  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1072  if (!track) continue;
1073  if(TMath::Abs(track->Eta())>0.9) continue;
1074  if (track->Pt()<0.15) continue;
1075  if (!(track->TestFilterBit(768))) continue;
1076 
1077  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1078  trackList->Add(track);
1079  triggers[iTT] = iTrack;
1080  iTT++;
1081 
1082  }
1083  }
1084  }
1085 
1086  if (iTT == 0) return -99999;
1087  Int_t nbRn = 0, index = 0 ;
1088  TRandom3* random = new TRandom3(0);
1089  nbRn = random->Integer(iTT);
1090 
1091  index = triggers[nbRn];
1092  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1093  return index;
1094 
1095 }
1096 
1097 //__________________________________________________________________________________
1099 
1100  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1101  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1102  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1103  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1104  double dphi = mphi-vphi;
1105  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1106  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1107  return dphi;//dphi in [-Pi, Pi]
1108 }
1109 
1110 
1111 //________________________________________________________________________
1113  //
1114  // retrieve event objects
1115  //
1117  return kFALSE;
1118 
1119  return kTRUE;
1120 }
1121 
1122 //_______________________________________________________________________
1124 {
1125  // Called once at the end of the analysis.
1126 
1127  AliInfo("Terminate");
1128  AliAnalysisTaskSE::Terminate();
1129 
1130  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1131  if (!fOutput) {
1132  AliError("fOutput not available");
1133  return;
1134  }
1135 
1136  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1137  if (!fTreeObservableTagging){
1138  Printf("ERROR: fTreeObservableTagging not available");
1139  return;
1140  }
1141 
1142 }
1143 
1144 //_________________________________________________________________________
1146 
1147  std::vector<fastjet::PseudoJet> fInputVectors;
1148  fInputVectors.clear();
1149  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1150  fastjet::PseudoJet PseudoTracks;
1151  fastjet::PseudoJet MyJet;
1152  fastjet::PseudoJet PseudoTracksCMS;
1153  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1154  //cout<<"CALL TO SOFTDROP"<<endl;
1155 
1156  Double_t zeta=0;
1157  Double_t angle=0;
1158  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1159  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1160  if (!fTrk) continue;
1161  JetInvMass += fTrk->M();
1162 
1163  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1164  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1165  TrackEnergy += fTrk->E();
1166  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1167  PseudJetInvMass += PseudoTracks.m();
1168  fInputVectors.push_back(PseudoTracks);
1169 
1170  }
1171 
1172  //fRandom->SetSeed(0);
1173  //here add N tracks with random phi and eta and theta according to bdmps distrib.
1174 
1175  MyJet.reset(fJet->Px(),fJet->Py(),fJet->Pz(),fJet->E());
1176  Double_t omegac=0.5*fqhat*fxlength*fxlength/0.2;
1177  Double_t thetac=TMath::Sqrt(12*0.2/(fqhat*TMath::Power(fxlength,3)));
1178  Double_t xQs=TMath::Sqrt(fqhat*fxlength);
1179  //cout<<"medium parameters "<<omegac<<" "<<thetac<<" "<<xQs<<endl;
1180 
1181  for(Int_t i=0;i<fAdditionalTracks;i++){
1182 
1183  Double_t ppx,ppy,ppz,kTscale,lim2o,lim1o;
1184  Double_t lim2=xQs;
1185  Double_t lim1=10000;
1186 
1187  //generation of kT according to 1/kT^4, with minimum QS=2 GeV and maximum ~sqrt(ptjet*T)
1188  fTf1Kt= new TF1("fTf1Kt","1/(x*x*x*x)",lim2,lim1);
1189  kTscale=fTf1Kt->GetRandom();
1190  //generation within the jet cone
1191 
1192  //generation of w according to 1/w, with minimum wc
1193  //omega needs to be larger than kT so to have well defined angles
1194  lim2o=kTscale;
1195  lim1o=kTscale/TMath::Sin(0.1);
1196  fTf1Omega= new TF1("fTf1Omega","1/x",lim2o,lim1o);
1197  Double_t omega=fTf1Omega->GetRandom();
1198 
1199  Double_t sinpptheta=kTscale/omega;
1200  Double_t pptheta=TMath::ASin(sinpptheta);
1201  //cout<<"angle_omega_kt"<<pptheta<<" "<<omega<<" "<<kTscale<<endl;
1202  if(pptheta>fJetRadius) continue;
1203 
1204  PseudoTracksCMS.reset(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1205  //boost the particle in the rest frame of the jet to the lab frame
1206  fastjet::PseudoJet PseudoTracksLab=PseudoTracksCMS.boost(MyJet);
1207  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1208  fInputVectors.push_back(PseudoTracksLab);
1209  //in the frame of the jet
1210  zeta=omega/fJet->E();
1211  angle=pptheta;
1212 
1213 
1214  }
1215 
1216 
1217 
1218 
1219 
1220 
1221 
1222 
1223 
1224  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1225  fastjet::contrib::Recluster *recluster;
1226  try {
1227  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1228  std::vector<fastjet::PseudoJet> fOutputJets;
1229  fOutputJets.clear();
1230  fOutputJets=fClustSeqSA.inclusive_jets(0);
1231 
1232  //cout<<fOutputJets[0].perp()<<" "<<fJet->Pt()<<endl;
1233 
1234  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1235 
1236  softdrop.set_verbose_structure(kTRUE);
1237  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1238  if(ReclusterAlgo==2) jetalgo=fastjet::antikt_algorithm;
1239  if(ReclusterAlgo==1) jetalgo=fastjet::kt_algorithm;
1240  if(ReclusterAlgo==0) jetalgo=fastjet::cambridge_algorithm;
1241 
1242  recluster = new fastjet::contrib::Recluster(jetalgo,1,true);
1243  softdrop.set_reclustering(true,recluster);
1244  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1245  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1246 
1247  Double_t SymParam, Mu, DeltaR, GroomedPt,GroomedMass;
1248  Int_t NGroomedBranches;
1249  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1250  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1251  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1252  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1253  GroomedPt=finaljet.perp();
1254  GroomedMass=finaljet.m();
1255  // if(beta==0){
1256  // if(ReclusterAlgo==0){
1257  // fShapesVar[12]=SymParam;
1258  // fShapesVar[13]=DeltaR;
1259  // fShapesVar[14]=zeta;
1260  // fShapesVar[15]=angle;
1261  // fShapesVar[16]=GroomedMass;}
1262  // if(ReclusterAlgo==1){
1263  // fShapesVar[17]=SymParam;
1264  // fShapesVar[18]=DeltaR;
1265  // fShapesVar[19]=zeta;
1266  // fShapesVar[20]=angle;
1267  // fShapesVar[21]=GroomedMass; }
1268 
1269  // if(ReclusterAlgo==2){
1270  // fShapesVar[22]=SymParam;
1271  // fShapesVar[23]=DeltaR;
1272  // fShapesVar[24]=zeta;
1273  // fShapesVar[25]=angle;
1274  // fShapesVar[26]=GroomedMass;
1275  // }}
1276  // if(beta==1){
1277  // fShapesVar[27]=SymParam;
1278  // fShapesVar[28]=DeltaR;
1279  // fShapesVar[29]=zeta;
1280  // fShapesVar[30]=angle;
1281  // fShapesVar[31]=GroomedMass;
1282  // }
1283  // //this one kills soft and large angle radiation
1284  // if((beta==1.5) && (zcut==0.5)){
1285  // fShapesVar[32]=SymParam;
1286  // fShapesVar[33]=DeltaR;
1287  // fShapesVar[34]=zeta;
1288  // fShapesVar[35]=angle;
1289  // fShapesVar[36]=GroomedMass; }
1290  // //this option favour democratic branches at large kt
1291  // if((beta==-1) && (zcut==0.005)){
1292  // fShapesVar[37]=SymParam;
1293  // fShapesVar[38]=DeltaR;
1294  // fShapesVar[39]=zeta;
1295  // fShapesVar[40]=angle;
1296  // fShapesVar[41]=GroomedMass; }
1297 
1298  // if((beta==-2) && (zcut==0.005)){
1299  // fShapesVar[42]=SymParam;
1300  // fShapesVar[43]=DeltaR;
1301  // fShapesVar[44]=zeta;
1302  // fShapesVar[45]=angle;
1303  // fShapesVar[46]=GroomedMass; }
1304 
1305 
1306 
1307 
1308 
1309 
1310 
1311 
1312 
1313  } catch (fastjet::Error) {
1314  AliError(" [w] FJ Exception caught.");
1315  //return -1;
1316  }
1317 
1318 
1319 
1320 
1321  if(recluster) { delete recluster; recluster = NULL; }
1322 
1323  if(fTf1Kt){ delete fTf1Kt;}
1324  if(fTf1Omega){ delete fTf1Omega;}
1325 
1326 
1327 
1328  return;
1329 
1330 
1331 }
1332 
1333 
1334 //_________________________________________________________________________
1336 
1337  std::vector<fastjet::PseudoJet> fInputVectors;
1338  fInputVectors.clear();
1339  fastjet::PseudoJet PseudoTracks;
1340  fastjet::PseudoJet PseudoTracksLab;
1341  double xflagalgo=0;
1342  double lnpt_relinject=0;
1343  double yinject=0;
1344  int xflagAdded=0;
1345  double zinject,angleinject,pptheta,sinpptheta,omega,omega2,angle2;
1346  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1347  Float_t pTscale=0., phiscale=0., thetascale=0., pXscale=0., pYscale=0., pZscale=0., pscale=0.;
1348 
1349  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1350  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1351  if (!fTrk) continue;
1352  if (fScaleELoss){
1353  pTscale = xfraction*sqrt(pow(fTrk->Px(),2)+pow(fTrk->Py(),2));
1354  phiscale = fTrk->Phi();
1355  thetascale = 2.*TMath::ATan(TMath::Exp(-1.*(fTrk->Eta())));
1356  pXscale = pTscale * TMath::Cos(phiscale);
1357  pYscale = pTscale * TMath::Sin(phiscale);
1358  pZscale = pTscale/TMath::Tan(thetascale);
1359  pscale = TMath::Sqrt(pTscale*pTscale+pZscale*pZscale);
1360  PseudoTracks.reset(pXscale, pYscale, pZscale, pscale);
1361  }
1362  else PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1363  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1364  fInputVectors.push_back(PseudoTracks);
1365 
1366  }
1367  if(fAddMedScat){
1368  for(int i = 0; i < fAddMedScatN; i++){
1369  TRandom3 rand1(0),rand2(0); //set range +- jet R
1370  Double_t randN1 = 0.4*0.4*rand1.Rndm();
1371  Double_t randN2 = 2*TMath::Pi()*rand2.Rndm();
1372  Double_t phi_rand = (fJet->Phi())+TMath::Sqrt(randN1)*TMath::Sin(randN2);
1373  Double_t eta_rand = (fJet->Eta())+TMath::Sqrt(randN1)*TMath::Cos(randN2);
1374  Double_t fAddMedScatPt = (fAddMedScatPtFrac*fJet->Pt())/fAddMedScatN;
1375  PseudoTracks.reset(fAddMedScatPt*TMath::Cos(phi_rand),fAddMedScatPt*TMath::Sin(phi_rand),fAddMedScatPt/TMath::Tan(eta_rand),fAddMedScatPt);
1376  fInputVectors.push_back(PseudoTracks);
1377  }
1378  }
1379 
1380 
1381  //add tracks to the jet prior to the reclusterer in case of iterative mapping of splittings
1382 
1383  Double_t omegac=0.5*fqhat*fxlength*fxlength/0.2;
1384  Double_t thetac=TMath::Sqrt(12*0.2/(fqhat*TMath::Power(fxlength,3)));
1385  Double_t xQs=TMath::Sqrt(fqhat*fxlength);
1386 
1387 
1388  for(Int_t i=0;i<fAdditionalTracks;i++){
1389 
1390  Double_t ppx,ppy,ppz,kTscale,lim2o,lim1o;
1391  Double_t lim2=xQs;
1392  Double_t lim1=10000;
1393 
1394  //generation of kT according to 1/kT^4, with minimum QS=2 GeV and maximum ~sqrt(ptjet*T)
1395  fTf1Kt= new TF1("fTf1Kt","1/(x*x*x*x)",lim2,lim1);
1396  kTscale=fTf1Kt->GetRandom();
1397  //generation within the jet cone
1398 
1399  //generation of w according to 1/w, with minimum wc
1400  //omega needs to be larger than kT so to have well defined angles
1401  lim2o=kTscale;
1402  lim1o=kTscale/TMath::Sin(0.1);
1403  fTf1Omega= new TF1("fTf1Omega","1/x",lim2o,lim1o);
1404  omega=fTf1Omega->GetRandom();
1405  sinpptheta=kTscale/omega;
1406  pptheta=TMath::ASin(sinpptheta);
1407  if(pptheta>fJetRadius) continue;
1408 
1409  //Lorentz vector in the frame where the jet moves along z axis
1410  TLorentzVector pTrackCMS(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1411  TVector3 MyJet(fJet->Px(),fJet->Py(),fJet->Pz());
1412  TVector3 direction = MyJet.Unit();
1413  //rotate the track to the jet frame
1414  pTrackCMS.RotateUz(direction);
1415 
1416  //add the rotated track to the jet
1417  PseudoTracksLab.reset(pTrackCMS.Px(),pTrackCMS.Py(),pTrackCMS.Pz(),pTrackCMS.E());
1418 
1419  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1420 
1421  omega2=PseudoTracksLab.perp();
1422  angle2=pTrackCMS.Angle(MyJet);
1423 
1424  fInputVectors.push_back(PseudoTracksLab);
1425  xflagAdded=1;
1426  }
1427 
1428 
1429 
1430 
1431  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1432  if(ReclusterAlgo==0){ xflagalgo=0.5;
1433  jetalgo=fastjet::kt_algorithm ;}
1434 
1435  if(ReclusterAlgo==1){ xflagalgo=1.5;
1436  jetalgo=fastjet::cambridge_algorithm;}
1437  if(ReclusterAlgo==2){ xflagalgo=2.5;
1438  jetalgo=fastjet::antikt_algorithm;}
1439 
1440  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1441 
1442  try {
1443  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1444  std::vector<fastjet::PseudoJet> fOutputJets;
1445  fOutputJets.clear();
1446  fOutputJets=fClustSeqSA.inclusive_jets(0);
1447 
1448  fastjet::PseudoJet jj;
1449  fastjet::PseudoJet j1;
1450  fastjet::PseudoJet j2;
1451  jj=fOutputJets[0];
1452  double ndepth=0;
1453  double nsd=0;
1454  double nall=0;
1455  fShapesVar[10]=fOutputJets[0].perp();
1456  while(jj.has_parents(j1,j2)){
1457 
1458  if(j1.perp() < j2.perp()) swap(j1,j2);
1459  double delta_R=j1.delta_R(j2);
1460  double z=j2.perp()/(j1.perp()+j2.perp());
1461  double y =log(1.0/delta_R);
1462  double lnpt_rel=log(z*delta_R);
1463  nall=nall+1;
1464  if(z>fHardCutoff){
1465  ndepth=ndepth+1;
1466  nsd=nsd+1;
1467  if(nsd==1 && ReclusterAlgo==1) fShapesVar[9]=z;
1468  Double_t LundEntries[6] = {y,lnpt_rel,fOutputJets[0].perp(),xflagalgo,partonFlavor,ndepth};
1469  fHLundIterative->Fill(LundEntries);}
1470  jj=j1;}
1471 
1472  if(ReclusterAlgo==1){
1473  fShapesVar[6]=nsd;
1474  fShapesVar[7]=nall;}
1475 
1476  if(fAdditionalTracks>0 && xflagAdded>0){
1477  zinject=omega2/fOutputJets[0].perp();
1478  angleinject=angle2;
1479 
1480  yinject =log(1.0/angleinject);
1481  lnpt_relinject=log(zinject*angleinject);
1482  Double_t LundEntriesInject[4] = {yinject,lnpt_relinject,fOutputJets[0].perp(),fJet->Pt()};
1483  fHLundIterativeInject->Fill(LundEntriesInject);}
1484 
1485 
1486  } catch (fastjet::Error) {
1487  AliError(" [w] FJ Exception caught.");
1488  //return -1;
1489  }
1490 
1491 
1492 
1493  if(fTf1Kt){ delete fTf1Kt;}
1494  if(fTf1Omega){ delete fTf1Omega;}
1495 
1496 
1497  return;
1498 
1499 
1500 }
1501 
1502 
1503 
1504 
void SetRadius(Double_t val)
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetFirstOrderSubtractedAngularity() const
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb=0)
double Double_t
Definition: External.C:58
void SetJetMinPt(Double_t val)
Double_t GetSecondOrderSubtractedSigma2() const
Definition: External.C:236
Float_t GetPartonEta6() const
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t FjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option, Double_t Beta_SD=0.0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Float_t GetPythiaEventWeight() const
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo, Float_t PartonFlavor)
Float_t GetPartonPhi7() const
Float_t CoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t E() const
Definition: AliEmcalJet.h:119
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb=0)
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
Double_t Px() const
Definition: AliEmcalJet.h:106
Float_t GetJetCoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
void SetRecombSheme(Int_t val)
void SetJetAlgorithm(Int_t val)
int Int_t
Definition: External.C:63
Double_t Eta() const
Double_t Pt() const
unsigned int UInt_t
Definition: External.C:33
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb=0)
float Float_t
Definition: External.C:68
Float_t fqhat
Random number generator.
Task to store and correlate the MC shapes.
TTree * fTreeObservableTagging
! Tree with tagging variables subtracted MC or true MC or raw
Double_t RelativePhi(Double_t mphi, Double_t vphi)
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.
Double_t GetSecondOrderSubtractedAngularity() const
Double_t GetSubjetFraction(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer)
Definition: Option.C:68
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Double_t fCent
!event centrality
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb=0)
void SetJetMaxEta(Double_t val)
AliEmcalJet * GetNextAcceptJet()
TF1 * fTf1Kt
to generate omega according to BDMPS tail
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 GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb=0)
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
AliEmcalJet * GetJet(Int_t index)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb=0)
Bool_t FillHistograms()
Function filling histograms.
Double_t GetFirstOrderSubtractedCircularity() const
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
AliEmcalList * fOutput
!output list
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
void NTValues(AliEmcalJet *jet, Int_t jetContNb, Float_t *nTFractions)
const Int_t nVar
Double_t GetSecondOrderSubtractedCircularity() const
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
const AliEmcalPythiaInfo * GetPythiaInfo() const
Double_t Pz() const
Definition: AliEmcalJet.h:108
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, Double_t zcut, Double_t beta, Int_t ReclusterAlgo)
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 Sigma2(AliEmcalJet *jet, Int_t jetContNb=0)
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
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.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const