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