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