AliPhysics  3337bb0 (3337bb0)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "AliEmcalJetFinder.h"
42 #include "AliAODEvent.h"
44 
45 using std::cout;
46 using std::endl;
47 
49 
50 //________________________________________________________________________
53  fContainer(0),
54  fMinFractionShared(0),
55  fJetShapeType(kGenShapes),
56  fJetShapeSub(kNoSub),
57  fJetSelection(kInclusive),
58  fPtThreshold(-9999.),
59  fRMatching(0.2),
60  fJetRadius(0.4),
61  fSubjetRadius(0.2),
62  fSelectedShapes(0),
63  fSwitchKtNSub(0),
64  fSwitchMinNSub(0),
65  fSwitchAktNSub(0),
66  fSwitchSDKtNSub(0),
67  fSwitchSDMinNSub(0),
68  fAdditionalTracks(0),
69  fminpTTrig(20.),
70  fmaxpTTrig(50.),
71  fangWindowRecoil(0.6),
72  fSemigoodCorrect(0),
73  fHolePos(0),
74  fHoleWidth(0),
75  fRandom(0),
76  fCentSelectOn(kTRUE),
77  fCentMin(0),
78  fCentMax(10),
79  fOneConstSelectOn(kFALSE),
80  fDerivSubtrOrder(0),
81  fPhiJetCorr6(0x0),
82  fPhiJetCorr7(0x0),
83  fEtaJetCorr6(0x0),
84  fEtaJetCorr7(0x0),
85  fPtJetCorr(0x0),
86  fPtJet(0x0),
87  fhpTjetpT(0x0),
88  fhPt(0x0),
89  fhPhi(0x0),
90  fNbOfConstvspT(0x0),
91  fTreeObservableTagging(0x0)
92 
93 {
94  for(Int_t i=0;i<33;i++){
95  fShapesVar[i]=0;}
96  SetMakeGeneralHistograms(kTRUE);
97 }
98 
99 //________________________________________________________________________
101  AliAnalysisTaskEmcalJet(name, kTRUE),
102  fContainer(0),
103  fMinFractionShared(0),
104  fJetShapeType(kGenShapes),
105  fJetShapeSub(kNoSub),
106  fJetSelection(kInclusive),
107  fPtThreshold(-9999.),
108  fRMatching(0.2),
109  fSelectedShapes(0),
110  fSwitchKtNSub(0),
111  fSwitchMinNSub(0),
112  fSwitchAktNSub(0),
113  fSwitchSDKtNSub(0),
114  fSwitchSDMinNSub(0),
115  fAdditionalTracks(0),
116  fminpTTrig(20.),
117  fmaxpTTrig(50.),
118  fangWindowRecoil(0.6),
119  fSemigoodCorrect(0),
120  fHolePos(0),
121  fHoleWidth(0),
122  fRandom(0),
123  fCentSelectOn(kTRUE),
124  fCentMin(0),
125  fCentMax(10),
126  fOneConstSelectOn(kFALSE),
127  fDerivSubtrOrder(0),
128  fPhiJetCorr6(0x0),
129  fPhiJetCorr7(0x0),
130  fEtaJetCorr6(0x0),
131  fEtaJetCorr7(0x0),
132  fPtJetCorr(0x0),
133  fPtJet(0x0),
134  fhpTjetpT(0x0),
135  fhPt(0x0),
136  fhPhi(0x0),
137  fNbOfConstvspT(0x0),
138  fTreeObservableTagging(0x0)
139 
140 {
141  // Standard constructor.
142 
143 
144  for(Int_t i=0;i<48;i++){
145  fShapesVar[i]=0;}
146 
148 
149  DefineOutput(1, TList::Class());
150  DefineOutput(2, TTree::Class());
151 
152 
153 }
154 
155 //________________________________________________________________________
157 {
159  delete fTreeObservableTagging;
161  }
162 
163  if(fRandom) delete fRandom;
164 }
165 
166 //________________________________________________________________________
168 {
169  // Create user output.
171 
172  Bool_t oldStatus = TH1::AddDirectoryStatus();
173  TH1::AddDirectory(oldStatus);
174 
175  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
176 
177  //TH1::AddDirectory(oldStatus);
178 
179  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
180  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
181 
182  const Int_t nVar = 54;
183 
184  TString *fShapesVarNames = new TString [nVar];
185 
186  fShapesVarNames[0] = "partonCode";
187  fShapesVarNames[1] = "ptJet";
188  fShapesVarNames[2] = "ptDJet";
189  fShapesVarNames[3] = "mJet";
190  fShapesVarNames[4] = "nbOfConst";
191  fShapesVarNames[5] = "angularity";
192  fShapesVarNames[6] = "Nsubjet1kt";
193  fShapesVarNames[7] = "Nsubjet2kt";
194  fShapesVarNames[8] = "Nsubjet1Min";
195  fShapesVarNames[9] = "Nsubjet2Min";
196  fShapesVarNames[10] = "DeltaRkt";
197  fShapesVarNames[11] = "DeltaRMin";
198  fShapesVarNames[12] = "SDSymm";
199  fShapesVarNames[13] = "SDDeltaR";
200  fShapesVarNames[14] = "SDGroomedFrac";
201  fShapesVarNames[15] = "SDGroomedN";
202  fShapesVarNames[16] = "SDMass";
203  fShapesVarNames[17] = "SDSymmkt";
204  fShapesVarNames[18] = "SDDeltaRkt";
205  fShapesVarNames[19] = "SDGroomedFrackt";
206  fShapesVarNames[20] = "SDGroomedNkt";
207  fShapesVarNames[21] = "SDMasskt";
208  fShapesVarNames[22] = "SDSymmAkt";
209  fShapesVarNames[23] = "SDDeltaRAkt";
210  fShapesVarNames[24] = "SDGroomedFracAkt";
211  fShapesVarNames[25] = "SDGroomedNAkt";
212  fShapesVarNames[26] = "SDMassAkt";
213  fShapesVarNames[27] = "SDSymmBeta1";
214  fShapesVarNames[28] = "SDDeltaRBeta1";
215  fShapesVarNames[29] = "SDGroomedFracBeta1";
216  fShapesVarNames[30] = "SDGroomedNBeta1";
217  fShapesVarNames[31] = "SDMassBeta1";
218  fShapesVarNames[32] = "SDSymmBeta15zcut05";
219  fShapesVarNames[33] = "SDDeltaRBeta15zcut05";
220  fShapesVarNames[34] = "SDGroomedFracBeta15zcut05";
221  fShapesVarNames[35] = "SDGroomedNBeta15zcut05";
222  fShapesVarNames[36] = "SDMassBeta15zcut05";
223  fShapesVarNames[37] = "SDSymmBetanegzcut025";
224  fShapesVarNames[38] = "SDDeltaRBetanegzcut025";
225  fShapesVarNames[39] = "SDGroomedFracBetanegzcut025";
226  fShapesVarNames[40] = "SDGroomedNBetanegzcut025";
227  fShapesVarNames[41] = "SDMassBetanegzcut025";
228  fShapesVarNames[42] = "SDSymmForm";
229  fShapesVarNames[43] = "SDDeltaRForm";
230  fShapesVarNames[44] = "SDGroomedFracForm";
231  fShapesVarNames[45] = "SDGroomedNForm";
232  fShapesVarNames[46] = "SDMassForm";
233  fShapesVarNames[47] = "weightPythia";
234  fShapesVarNames[48] = "SDNsubjet1kt";
235  fShapesVarNames[49] = "SDNsubjet2kt";
236  fShapesVarNames[50] = "SDNsubjet1Min";
237  fShapesVarNames[51] = "SDNsubjet2Min";
238  fShapesVarNames[52] = "SDDeltaRkt_Nsub";
239  fShapesVarNames[53] = "SDDeltaRMin_Nsub";
240 
241 
242 
243  //fShapesVarNames[7] = "lesub";
244  //fShapesVarNames[8] = "CoreFraction";
245  //fShapesVarNames[9] = "Nsubjet1";
246  //fShapesVarNames[10] = "Nsubjet2";
247  //fShapesVarNames[11] = "DeltaR";
248  //fShapesVarNames[12] = "OpenAngle";
249  //fShapesVarNames[13] = "weightPythia";
250 
251  //fShapesVarNames[14] = "NT70";
252  //fShapesVarNames[15] = "nConstNT70";
253  //fShapesVarNames[16] = "NT80";
254  //fShapesVarNames[17] = "nConstNT80";
255  //fShapesVarNames[18] = "NT90";
256  //fShapesVarNames[19] = "nConstNT90";
257  //fShapesVarNames[20] = "NT95";
258  //fShapesVarNames[21] = "nConstNT95";
259 
260  //fShapesVarNames[22] = "SubjetFraction";
261 
262 
263  for(Int_t ivar=0; ivar < nVar; ivar++){
264  cout<<"looping over variables"<<endl;
265  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
266 
267 
268 
269  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
270  fOutput->Add(fPhiJetCorr6);
271  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
272  fOutput->Add(fEtaJetCorr6);
273 
274  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
275  fOutput->Add(fPhiJetCorr7);
276  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
277  fOutput->Add(fEtaJetCorr7);
278 
279  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
280  fOutput->Add(fPtJetCorr);
281  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
282  fOutput->Add(fPtJet);
283 
284  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
285  fOutput->Add(fhpTjetpT);
286  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
287  fOutput->Add(fhPt);
288  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
289  fOutput->Add(fhPhi);
290 
291  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
292  fOutput->Add(fNbOfConstvspT);
293 
294  //fOutput->Add(fTreeObservableTagging);
295 
296  fRandom = new TRandom3(0);
297  PostData(1, fOutput); // Post data for ALL output slots > 0 here
298  PostData(2, fTreeObservableTagging);
299 
300  delete [] fShapesVarNames;
301 
302  }
303 
304 //________________________________________________________________________
306 {
307  // Run analysis code here, if needed. It will be executed before FillHistograms().
308  // if (gRandom) delete gRandom;
309  // gRandom = new TRandom3(0);
310  return kTRUE;
311 }
312 
313 //________________________________________________________________________
315 {
316  // Fill histograms.
317  //cout<<"IntoFillHistograms"<<endl;
318  AliEmcalJet* jet1 = NULL;
319  AliJetContainer *jetCont = GetJetContainer(0);
320 
321  Float_t kWeight=1;
322  if (fCentSelectOn)
323  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
324 
325  AliAODTrack *triggerHadron = 0x0;
326 
327  if (fJetSelection == kRecoil) {
328  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
329  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
330 
331 
332  if (triggerHadronLabel==-99999) {
333  //Printf ("Trigger Hadron not found, return");
334  return 0;}
335 
336 
338  TClonesArray *trackArrayAn = partContAn->GetArray();
339  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
340 
341  if (!triggerHadron) {
342  //Printf("No Trigger hadron with the found label!!");
343  return 0;
344  }
345 
346  if(fSemigoodCorrect){
347  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
348  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
349  return 0;}
350  }
351 
352  fhPt->Fill(triggerHadron->Pt());
353 
354  }
355 
356  if(jetCont) {
357  jetCont->ResetCurrentID();
358  while((jet1 = jetCont->GetNextAcceptJet())) {
359  //Printf("jet1=%p", jet1);
360  if (!jet1) continue;
361  AliEmcalJet* jet2 = 0x0;
362  AliEmcalJet* jet3 = 0x0;
363  fPtJet->Fill(jet1->Pt());
364  AliEmcalJet *jetUS = NULL;
365  Int_t ifound=0, jfound=0;
366  Int_t ilab=-1, jlab=-1;
367 
369  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
370  if(TMath::Abs(disthole)<fHoleWidth){
371  continue;
372  }
373  }
374 
375  Float_t dphiRecoil = 0.;
376  if (fJetSelection == kRecoil){
377  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
378  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
379  // Printf("Recoil jets back to back not found! continuing");
380  continue;
381  }
382 
383  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
384  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
385  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
386 
387  }
388 
389 
390  fShapesVar[0] = 0.;
391 
392  if (fJetShapeType == kGenShapes){
393  const AliEmcalPythiaInfo *partonsInfo = 0x0;
394  partonsInfo = GetPythiaInfo();
395  //Printf("partonsInfo=%p", partonsInfo);
396  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
397  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
398  kWeight=partonsInfo->GetPythiaEventWeight();
399  //Printf("kWeight=%f", kWeight);
400  fShapesVar[47] = kWeight;
401 
402  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
403  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
404  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
405  if(dRp1 < fRMatching) {
406  fShapesVar[0] = partonsInfo->GetPartonFlag6();
407  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
408  }
409  else {
410  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
411  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
412  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
413  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
414  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
415  if(dRp1 < fRMatching) {
416  fShapesVar[0] = partonsInfo->GetPartonFlag7();
417  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
418  }
419  else fShapesVar[0]=0;
420  }
421  }
422 
423  Double_t ptSubtracted = 0;
424  ptSubtracted= jet1->Pt();
425  //Printf("ptSubtracted=%f", ptSubtracted);
426 
427 
428  if (ptSubtracted < fPtThreshold) continue;
429 
430  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() < 2)) continue;
431 
432  //AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
433  //Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
434 
435 
436 
437 
438 
439 
440 
441  fShapesVar[1] = ptSubtracted;
442  fShapesVar[2] = GetJetpTD(jet1,0);
443  fShapesVar[3] = GetJetMass(jet1,0);
444  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
445  fShapesVar[5] = GetJetAngularity(jet1,0);
446  //nsub1 and nsub2 for kT
447  fShapesVar[6] = 0;
448  if(fSwitchKtNSub==1) fShapesVar[6]=FjNSubJettiness(jet1,0,1,0,1,0);
449  fShapesVar[7] = 0;
450  if(fSwitchKtNSub==1) fShapesVar[7]=FjNSubJettiness(jet1,0,2,0,1,0);
451 
452  //nsub1 and nsub2 for min_axis
453  fShapesVar[8] =0;
454  if(fSwitchMinNSub==1) fShapesVar[8]=FjNSubJettiness(jet1,0,1,10,1,0);
455  fShapesVar[9] = 0;
456  if(fSwitchMinNSub==1) fShapesVar[9]=FjNSubJettiness(jet1,0,2,10,1,0);
457  //nsub1 and nsub2 for akt
458  fShapesVar[10] = 0;
459  if(fSwitchKtNSub==1) fShapesVar[10]=FjNSubJettiness(jet1,0,2,0,1,1);
460  fShapesVar[11] =0;
461  if(fSwitchMinNSub==1) fShapesVar[11]=FjNSubJettiness(jet1,0,2,10,1,1);
462  //nsub1 and nsub2 for kt with SD with Beta = 0 and Zcut =0.1
463  fShapesVar[48] =0;
464  if(fSwitchSDKtNSub==1) fShapesVar[48]=FjNSubJettiness(jet1,0,1,0,1,0,0,0.1,1);
465  fShapesVar[49] =0;
466  if(fSwitchSDKtNSub==1) fShapesVar[49]=FjNSubJettiness(jet1,0,2,0,1,0,0,0.1,1);
467  //nsub1 and nsub2 for min_axis with SD with Beta = 0 and Zcut =0.1
468  fShapesVar[50] =0;
469  if(fSwitchSDMinNSub==1) fShapesVar[50]=FjNSubJettiness(jet1,0,1,10,1,0,0,0.1,1);
470  fShapesVar[51] =0;
471  if(fSwitchSDMinNSub==1) fShapesVar[51]=FjNSubJettiness(jet1,0,2,10,1,0,0,0.1,1);
472  //deltaR using axes from 2 subjettiness with kt and Min algorithms with soft drop
473  fShapesVar[52] =0;
474  if(fSwitchSDKtNSub==1) fShapesVar[52]=FjNSubJettiness(jet1,0,2,0,1,1,0,0.1,1);
475  fShapesVar[53] =0;
476  if(fSwitchSDMinNSub==1) fShapesVar[53]=FjNSubJettiness(jet1,0,2,10,1,1,0,0.1,1);
477 
478  //SoftDropParameters for different reclustering strategies and beta values
479  SoftDrop(jet1,jetCont,0.1,0,0);
480  SoftDrop(jet1,jetCont,0.1,0,1);
481  SoftDrop(jet1,jetCont,0.1,0,2);
482  SoftDrop(jet1,jetCont,0.1,1,0);
483  SoftDrop(jet1,jetCont,0.5,1.5,0);
484  SoftDrop(jet1,jetCont,0.25,-1,0);
485  SoftDrop(jet1,jetCont,0.1,-2,0);
486  // Float_t nTFractions[8]={0.,0.,0.,0.,0.,0.,0.,0.};
487  //NTValues(jet1, 0, nTFractions);
488  //shape 13 is pythia weight!
489  //for (Int_t ishape=14; ishape<22; ishape++) fShapesVar[ishape] = nTFractions[ishape-14];
490 
491  //fShapesVar[22]= GetSubjetFraction(jet1,0,fJetRadius,Reclusterer1);
492 
493  fTreeObservableTagging->Fill();
494 
495  }
496 
497  }
498 
499  return kTRUE;
500 }
501 
502 //________________________________________________________________________
504  //calc subtracted jet mass
505  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
507  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
508  else
509  return jet->M();
510 }
511 
512 //________________________________________________________________________
514 
515  AliJetContainer *jetCont = GetJetContainer(jetContNb);
516  if (!jet->GetNumberOfTracks())
517  return 0;
518  Double_t den=0.;
519  Double_t num = 0.;
520  AliVParticle *vp1 = 0x0;
521  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
522  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
523 
524  if (!vp1){
525  Printf("AliVParticle associated to constituent not found");
526  continue;
527  }
528 
529  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
530  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
531  Double_t dr = TMath::Sqrt(dr2);
532  num=num+vp1->Pt()*dr;
533  den=den+vp1->Pt();
534  }
535  return num/den;
536 }
537 
538 //________________________________________________________________________
540 
541  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
544  else
545  return Angularity(jet, jetContNb);
546 
547 }
548 
549 
550 //________________________________________________________________________
552 
553  AliJetContainer *jetCont = GetJetContainer(jetContNb);
554  if (!jet->GetNumberOfTracks())
555  return 0;
556  Double_t den=0.;
557  Double_t num = 0.;
558  AliVParticle *vp1 = 0x0;
559  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
560  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
561 
562  if (!vp1){
563  Printf("AliVParticle associated to constituent not found");
564  continue;
565  }
566 
567  num=num+vp1->Pt()*vp1->Pt();
568  den=den+vp1->Pt();
569  }
570  return TMath::Sqrt(num)/den;
571 }
572 
573 //________________________________________________________________________
575  //calc subtracted jet mass
576  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
578  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
579  else
580  return PTD(jet, jetContNb);
581 
582 }
583 
584 //_____________________________________________________________________________
586 
587  AliJetContainer *jetCont = GetJetContainer(jetContNb);
588  if (!jet->GetNumberOfTracks())
589  return 0;
590  Double_t mxx = 0.;
591  Double_t myy = 0.;
592  Double_t mxy = 0.;
593  int nc = 0;
594  Double_t sump2 = 0.;
595  Double_t pxjet=jet->Px();
596  Double_t pyjet=jet->Py();
597  Double_t pzjet=jet->Pz();
598 
599 
600  //2 general normalized vectors perpendicular to the jet
601  TVector3 ppJ1(pxjet, pyjet, pzjet);
602  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
603  ppJ3.SetMag(1.);
604  TVector3 ppJ2(-pyjet, pxjet, 0);
605  ppJ2.SetMag(1.);
606  AliVParticle *vp1 = 0x0;
607  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
608  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
609 
610  if (!vp1){
611  Printf("AliVParticle associated to constituent not found");
612  continue;
613  }
614 
615  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
616 
617  //local frame
618  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
619  TVector3 pPerp = pp - pLong;
620  //projection onto the two perpendicular vectors defined above
621 
622  Float_t ppjX = pPerp.Dot(ppJ2);
623  Float_t ppjY = pPerp.Dot(ppJ3);
624  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
625  if(ppjT<=0) return 0;
626 
627  mxx += (ppjX * ppjX / ppjT);
628  myy += (ppjY * ppjY / ppjT);
629  mxy += (ppjX * ppjY / ppjT);
630  nc++;
631  sump2 += ppjT;}
632 
633  if(nc<2) return 0;
634  if(sump2==0) return 0;
635  // Sphericity Matrix
636  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
637  TMatrixDSym m0(2,ele);
638 
639  // Find eigenvectors
640  TMatrixDSymEigen m(m0);
641  TVectorD eval(2);
642  TMatrixD evecm = m.GetEigenVectors();
643  eval = m.GetEigenValues();
644  // Largest eigenvector
645  int jev = 0;
646  // cout<<eval[0]<<" "<<eval[1]<<endl;
647  if (eval[0] < eval[1]) jev = 1;
648  TVectorD evec0(2);
649  // Principle axis
650  evec0 = TMatrixDColumn(evecm, jev);
651  Double_t compx=evec0[0];
652  Double_t compy=evec0[1];
653  TVector2 evec(compx, compy);
654  Double_t circ=0;
655  if(jev==1) circ=2*eval[0];
656  if(jev==0) circ=2*eval[1];
657 
658  return circ;
659 
660 
661 
662 }
663 
664 
665 
666 
667 //________________________________________________________________________
669  //calc subtracted jet mass
670 
671  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
674  else
675  return Circularity(jet, jetContNb);
676 
677 }
678 
679 //________________________________________________________________________
681 
682  AliJetContainer *jetCont = GetJetContainer(jetContNb);
683  if (!jet->GetNumberOfTracks())
684  return 0;
685  Double_t den=0.;
686  Double_t num = 0.;
687  AliVParticle *vp1 = 0x0;
688  AliVParticle *vp2 = 0x0;
689  std::vector<int> ordindex;
690  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
691  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
692  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
693 
694  if(ordindex.size()<2) return -1;
695 
696  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
697  if (!vp1){
698  Printf("AliVParticle associated to Leading constituent not found");
699  return -1;
700  }
701 
702  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
703  if (!vp2){
704  Printf("AliVParticle associated to Subleading constituent not found");
705  return -1;
706  }
707 
708 
709  num=vp1->Pt();
710  den=vp2->Pt();
711  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
712 
713 return num-den;
714 }
715 
716 //________________________________________________________________________
718  //calc subtracted jet mass
719 
720  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
723  else
724  return LeSub(jet, jetContNb);
725 
726 }
727 
728 //________________________________________________________________________
730  //calc subtracted jet mass
731 
732  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
735  else
736  return jet->GetNumberOfTracks();
737 
738  }
739 
740 
741 //________________________________________________________________________
742 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
743 
744  AliJetContainer *JetCont = GetJetContainer(JetContNb);
745  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
746  Reclusterer->SetRadius(SubJetRadius);
747  Reclusterer->SetJetMinPt(SubJetMinPt);
748  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
749  Reclusterer->SetJetMaxEta(0.9-JetRadius);
750  Reclusterer->SetRecombSheme(0);
751  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
752 
753  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
754  Double_t dVtx[3]={0.,0.,0.};
755  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
756  return Reclusterer;
757 }
758 
759 
760 
761 
762 //----------------------------------------------------------------------
764  AliJetContainer *JetCont = GetJetContainer(JetContNb);
765  AliEmcalJet *SubJet=NULL;
766  Double_t DeltaR1=0;
767  Double_t DeltaR2=0;
768  AliVParticle *JetParticle=0x0;
769  Double_t SubJetiness_Numerator = 0;
770  Double_t SubJetiness_Denominator = 0;
771  Double_t Index=-2;
772  if (Reclusterer->GetNumberOfJets() < 1) return -2;
773  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
774  if(Index==-999) return -2;
775  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
776  SubJetiness_Denominator=Jet->Pt();
777  return SubJetiness_Numerator/SubJetiness_Denominator;
778 
779 
780 }
781 //__________________________________________________________________________________
783 
784  AliJetContainer *jetCont = GetJetContainer(jetContNb);
785  if (!jet->GetNumberOfTracks())
786  return 0;
787  Double_t den=0.;
788  Double_t num = 0.;
789  AliVParticle *vp1 = 0x0;
790  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
791  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
792 
793  if (!vp1){
794  Printf("AliVParticle associated to constituent not found");
795  continue;
796  }
797 
798  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
799  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
800  Double_t dr = TMath::Sqrt(dr2);
801  if(dr<=fSubjetRadius) num=num+vp1->Pt();
802 
803  }
804  return num/jet->Pt();
805 }
806 
807 
808 
809 
810 //________________________________________________________________________
812  //calc subtracted jet mass
813 
814  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
817  else
818  return CoreFrac(jet, jetContNb);
819 
820 }
821 
822 
823 
824 
825 //----------------------------------------------------------------------
827  AliEmcalJet *SubJet=NULL;
828  Double_t SortingVariable;
829  Int_t ArraySize =N+1;
830  TArrayD *JetSorter = new TArrayD(ArraySize);
831  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
832  for (Int_t i=0; i<ArraySize; i++){
833  JetSorter->SetAt(0,i);
834  }
835  for (Int_t i=0; i<ArraySize; i++){
836  JetIndexSorter->SetAt(0,i);
837  }
838  if(Reclusterer->GetNumberOfJets()<N) return -999;
839  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
840  SubJet=Reclusterer->GetJet(i);
841  if (Type==0) SortingVariable=SubJet->Pt();
842  else if (Type==1) SortingVariable=SubJet->E();
843  else if (Type==2) SortingVariable=SubJet->M();
844  for (Int_t j=0; j<N; j++){
845  if (SortingVariable>JetSorter->GetAt(j)){
846  for (Int_t k=N-1; k>=j; k--){
847  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
848  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
849  }
850  JetSorter->SetAt(SortingVariable,j);
851  JetIndexSorter->SetAt(i,j);
852  break;
853  }
854  }
855  }
856  if (!Index) return JetSorter->GetAt(N-1);
857  else return JetIndexSorter->GetAt(N-1);
858 }
859 
860 
861 
862 //returns -1 if the Nth hardest jet is requested where N>number of available jets
863 //type: 0=Pt 1=E 2=M
864 //Index TRUE=returns index FALSE=returns value of quantatiy in question
865 
866 
867 
868 
869 
870 //______________________________________________________________________________
872 
873  AliJetContainer *jetCont = GetJetContainer(jetContNb);
874  if (!jet->GetNumberOfTracks())
875  return 0;
876  Double_t mxx = 0.;
877  Double_t myy = 0.;
878  Double_t mxy = 0.;
879  int nc = 0;
880  Double_t sump2 = 0.;
881 
882  AliVParticle *vp1 = 0x0;
883  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
884  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
885 
886  if (!vp1){
887  Printf("AliVParticle associated to constituent not found");
888  continue;
889  }
890 
891  Double_t ppt=vp1->Pt();
892  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
893 
894  Double_t deta = vp1->Eta()-jet->Eta();
895  mxx += ppt*ppt*deta*deta;
896  myy += ppt*ppt*dphi*dphi;
897  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
898  nc++;
899  sump2 += ppt*ppt;
900 
901  }
902  if(nc<2) return 0;
903  if(sump2==0) return 0;
904  // Sphericity Matrix
905  Double_t ele[4] = {mxx , mxy , mxy , myy };
906  TMatrixDSym m0(2,ele);
907 
908  // Find eigenvectors
909  TMatrixDSymEigen m(m0);
910  TVectorD eval(2);
911  TMatrixD evecm = m.GetEigenVectors();
912  eval = m.GetEigenValues();
913  // Largest eigenvector
914  int jev = 0;
915  // cout<<eval[0]<<" "<<eval[1]<<endl;
916  if (eval[0] < eval[1]) jev = 1;
917  TVectorD evec0(2);
918  // Principle axis
919  evec0 = TMatrixDColumn(evecm, jev);
920  Double_t compx=evec0[0];
921  Double_t compy=evec0[1];
922  TVector2 evec(compx, compy);
923  Double_t sig=0;
924  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
925  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
926 
927  return sig;
928 
929 }
930 
931 //________________________________________________________________________
933  //calc subtracted jet mass
934 
935  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
938  else
939  return Sigma2(jet, jetContNb);
940 
941 }
942 
943 
944 //_________________________________________________________________________
946 
947  AliJetContainer *jetCont = GetJetContainer(jetContNb);
948  if (!jet->GetNumberOfTracks())
949  return;
950 
951  Double_t ptJet = jet->Pt();
952 
953  AliVParticle *vp1 = 0x0;
954  std::vector<int> ordindex;
955  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
956  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
957  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
958  //if(ordindex.size()<2) return -1;
959 
960  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
961 
962  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
963  if (!vp1){
964  Printf("AliVParticle associated to Leading constituent not found");
965  return;
966  }
967 
968  if (nTFractions[0] <= 0.7*ptJet){
969  nTFractions[0] += vp1->Pt();
970  nTFractions[1] +=1;
971  }
972 
973  if (nTFractions[2] <= 0.8*ptJet){
974  nTFractions[2] += vp1->Pt();
975  nTFractions[3] +=1;
976  }
977 
978  if (nTFractions[4] <= 0.9*ptJet){
979  nTFractions[4] += vp1->Pt();
980  nTFractions[5] +=1;
981  }
982 
983  if (nTFractions[6] <= 0.95*ptJet){
984  nTFractions[6] += vp1->Pt();
985  nTFractions[7] +=1;
986  }
987  }
988 }
989 //_________________________________________________________________________________________________
991 
992  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
993 
994  //Algorithm==0 -> kt_axes;
995  // Algorithm==1 -> ca_axes;
996  //Algorithm==2 -> antikt_0p2_axes;
997  //Algorithm==3 -> wta_kt_axes;
998  //Algorithm==4 -> wta_ca_axes;
999  //Algorithm==5 -> onepass_kt_axes;
1000  //Algorithm==6 -> onepass_ca_axes;
1001  //Algorithm==7 -> onepass_antikt_0p2_axes;
1002  //Algorithm==8 -> onepass_wta_kt_axes;
1003  //Algorithm==9 -> onepass_wta_ca_axes;
1004  //Algorithm==10 -> min_axes;
1005 
1006 
1007  //Option==0 returns Nsubjettiness Value
1008  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1009  //Option==2 && N==2 returns Delta R
1010 
1011  if (Jet->GetNumberOfTracks()>=N){
1012  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1013  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1014  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1015  JetFinder->SetRadius(fJetRadius);
1016  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1017  JetFinder->SetRecombSheme(0);
1018  JetFinder->SetJetMinPt(Jet->Pt());
1019  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1020  Double_t dVtx[3]={1,1,1};
1021  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
1022  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,0.2,Beta,Option,0,Beta_SD,ZCut,SoftDropOn);
1023 
1024  }
1025  else return -2;
1026 }
1027 
1028 
1029 
1030 //________________________________________________________________________
1032 
1034  TClonesArray *tracksArray = partCont->GetArray();
1035 
1036  if(!partCont || !tracksArray) return -99999;
1037  AliAODTrack *track = 0x0;
1038  AliEmcalParticle *emcPart = 0x0;
1039 
1040 
1041  TList *trackList = new TList();
1042  Int_t triggers[100];
1043  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1044  Int_t iTT = 0;
1045 
1046  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
1047 
1048 
1049  if (fJetShapeSub == kConstSub){
1050  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
1051  if (!emcPart) continue;
1052  if(TMath::Abs(emcPart->Eta())>0.9) continue;
1053  if (emcPart->Pt()<0.15) continue;
1054 
1055  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
1056  trackList->Add(emcPart);
1057  triggers[iTT] = iTrack;
1058  iTT++;
1059  }
1060  }
1061  else{
1062  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1063  if (!track) continue;
1064  if(TMath::Abs(track->Eta())>0.9) continue;
1065  if (track->Pt()<0.15) continue;
1066  if (!(track->TestFilterBit(768))) continue;
1067 
1068  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1069  trackList->Add(track);
1070  triggers[iTT] = iTrack;
1071  iTT++;
1072 
1073  }
1074  }
1075  }
1076 
1077  if (iTT == 0) return -99999;
1078  Int_t nbRn = 0, index = 0 ;
1079  TRandom3* random = new TRandom3(0);
1080  nbRn = random->Integer(iTT);
1081 
1082  index = triggers[nbRn];
1083  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1084  return index;
1085 
1086 }
1087 
1088 //__________________________________________________________________________________
1090 
1091  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1092  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1093  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1094  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1095  double dphi = mphi-vphi;
1096  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1097  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1098  return dphi;//dphi in [-Pi, Pi]
1099 }
1100 
1101 
1102 //________________________________________________________________________
1104  //
1105  // retrieve event objects
1106  //
1108  return kFALSE;
1109 
1110  return kTRUE;
1111 }
1112 
1113 //_______________________________________________________________________
1115 {
1116  // Called once at the end of the analysis.
1117 
1118  AliInfo("Terminate");
1119  AliAnalysisTaskSE::Terminate();
1120 
1121  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1122  if (!fOutput) {
1123  AliError("fOutput not available");
1124  return;
1125  }
1126 
1127  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1128  if (!fTreeObservableTagging){
1129  Printf("ERROR: fTreeObservableTagging not available");
1130  return;
1131  }
1132 
1133 }
1134 
1135 //_________________________________________________________________________
1136 void AliAnalysisTaskEmcalJetShapesMC::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, int ReclusterAlgo){
1137 
1138  std::vector<fastjet::PseudoJet> fInputVectors;
1139  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1140 
1141  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1142  //cout<<"CALL TO SOFTDROP"<<endl;
1143  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
1144  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
1145  UShort_t FJNTracks=0,EmcalJetNTracks=0;
1146 
1147  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1148  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1149  if (!fTrk) continue;
1150  JetInvMass += fTrk->M();
1151 
1152  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1153  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1154  TrackEnergy += fTrk->E();
1155  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1156  PseudJetInvMass += PseudoTracks.m();
1157  fInputVectors.push_back(PseudoTracks);
1158  EmcalJetTrackEta[i]=fTrk->Eta();
1159  EmcalJetTrackPhi[i]=fTrk->Phi();
1160  EmcalJetTrackPt[i]=fTrk->Pt();
1161  EmcalJetNTracks++;
1162  }
1163 
1164  fRandom->SetSeed(0);
1165  //here add N tracks with random phi and eta and theta according to bdmps distrib.
1166  fastjet::PseudoJet MyJet(fJet->Px(),fJet->Py(),fJet->Pz(),fJet->E());
1167 
1168  for(Int_t i=0;i<fAdditionalTracks;i++){
1169  Double_t ppx,ppy,ppz,kTscale,pTscale,ppE;
1170  Double_t lim2=2;
1171  Double_t lim1=10;
1172  Double_t xpower=-4;
1173  Double_t xpowero=-2;
1174 
1175 
1176  Double_t part1=0;
1177  Double_t part2=0;
1178  Double_t part1o=0;
1179  Double_t part2o=0;
1180 
1181  //coordinates in the reference frame of the jet
1182  Double_t ppphi =2*TMath::Pi()*fRandom->Uniform();
1183  //generation of kT according to 1/kT^4, with minimum QS=2 GeV
1184  part1=(TMath::Power(lim1,xpower+1)-TMath::Power(lim2,xpower+1))*fRandom->Uniform();
1185  part2=TMath::Power(lim2,xpower+1);
1186  kTscale=TMath::Power(part1+part2,1/(xpower+1));
1187 
1188 
1189 
1190 
1191  //generation of w according to 1/w2, with minimum wc
1192  //omega needs to be larger than kT so to have well defined angles
1193  Double_t lim1o=10;
1194  Double_t lim2o=kTscale;
1195  part1o=(TMath::Power(lim1o,xpowero+1)-TMath::Power(lim2o,xpowero+1))*fRandom->Uniform();
1196  part2o=TMath::Power(lim2o,xpowero+1);
1197  Double_t omega=TMath::Power(part1o+part2o,1/(xpowero+1));
1198 
1199  Double_t sinpptheta=kTscale/omega;
1200  Double_t pptheta=TMath::ASin(sinpptheta);
1201 
1202  if(pptheta>0.4) continue;
1203  fastjet::PseudoJet PseudoTracksCMS(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1204  //cout<<PseudoTracksCMS.px()<<" "<<PseudoTracksCMS.py()<<" "<<PseudoTracksCMS.perp()<<" "<<PseudoTracksCMS.e()<<endl;
1205  //boost the particle to the lab frame
1206  fastjet::PseudoJet PseudoTracksLab=PseudoTracksCMS.boost(MyJet);
1207  //cout<<" "<<PseudoTracksLab.px()<<" "<<PseudoTracksLab.py()<<" "<<PseudoTracksLab.pz()<<" "<<PseudoTracksLab.e()<<" "<<PseudoTracksLab.perp()<<endl;
1208  //cout<<" "<<fJet->Px()<<" "<<fJet->Py()<<" "<<fJet->Pz()<<" "<<fJet->E()<<" "<<fJet->Pt()<<endl;
1209 
1210  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1211  fInputVectors.push_back(PseudoTracksLab);}
1212 
1213 
1214 
1215 
1216  fastjet::JetDefinition *fJetDef;
1217  fastjet::ClusterSequence *fClustSeqSA;
1218 
1219 
1220 
1221  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1222 
1223  try {
1224  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
1225  } catch (fastjet::Error) {
1226  AliError(" [w] FJ Exception caught.");
1227  //return -1;
1228  }
1229 
1230  std::vector<fastjet::PseudoJet> fOutputJets;
1231  fOutputJets.clear();
1232  fOutputJets=fClustSeqSA->inclusive_jets(0);
1233 
1234 
1235  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
1236  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1237  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
1238  softdrop.set_verbose_structure(kTRUE);
1239  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
1240  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
1241  fastjet::contrib::Recluster *recluster;
1242  if(ReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
1243  if(ReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
1244  if(ReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
1245  softdrop.set_reclustering(true,recluster);
1246  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1247  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
1248  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1249  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1250 
1251 
1252  //AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
1253  //std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
1254  //Double_t FastjetTrackDelR,EmcalTrackDelR;
1255  //for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
1256  // if(i<=finaljet.constituents().size()){
1257  // FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
1258  // FJTrackEta[i]=fSDTracks[i].eta();
1259  // FJTrackPhi[i]=fSDTracks[i].phi();
1260  // FJTrackPt[i]=fSDTracks[i].perp();
1261  // FJNTracks++;
1262  // }
1263  // AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1264  // EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
1265  //}
1266  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1267  //Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
1268  //jet->SetNumberOfTracks(nConstituents);
1269  Double_t SymParam, Mu, DeltaR, GroomedPt,GroomedMass;
1270  Int_t NGroomedBranches;
1271  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1272  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1273  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1274  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1275  GroomedPt=finaljet.perp();
1276  GroomedMass=finaljet.m();
1277  if(beta==0){
1278  if(ReclusterAlgo==0){
1279  fShapesVar[12]=SymParam;
1280  fShapesVar[13]=DeltaR;
1281  fShapesVar[14]=GroomedPt;
1282  fShapesVar[15]=NGroomedBranches;
1283  fShapesVar[16]=GroomedMass;}
1284  if(ReclusterAlgo==1){
1285  fShapesVar[17]=SymParam;
1286  fShapesVar[18]=DeltaR;
1287  fShapesVar[19]=GroomedPt;
1288  fShapesVar[20]=NGroomedBranches;
1289  fShapesVar[21]=GroomedMass; }
1290 
1291  if(ReclusterAlgo==2){
1292  fShapesVar[22]=SymParam;
1293  fShapesVar[23]=DeltaR;
1294  fShapesVar[24]=GroomedPt;
1295  fShapesVar[25]=NGroomedBranches;
1296  fShapesVar[26]=GroomedMass;
1297  }}
1298  if(beta==1){
1299  fShapesVar[27]=SymParam;
1300  fShapesVar[28]=DeltaR;
1301  fShapesVar[29]=GroomedPt;
1302  fShapesVar[30]=NGroomedBranches;
1303  fShapesVar[31]=GroomedMass;
1304  }
1305  //this one kills soft and large angle radiation
1306  if((beta==1.5) && (zcut==0.5)){
1307  fShapesVar[32]=SymParam;
1308  fShapesVar[33]=DeltaR;
1309  fShapesVar[34]=GroomedPt;
1310  fShapesVar[35]=NGroomedBranches;
1311  fShapesVar[36]=GroomedMass; }
1312  //this option favour democratic branches at large kt
1313  if((beta==-1) && (zcut==0.25)){
1314  fShapesVar[37]=SymParam;
1315  fShapesVar[38]=DeltaR;
1316  fShapesVar[39]=GroomedPt;
1317  fShapesVar[40]=NGroomedBranches;
1318  fShapesVar[41]=GroomedMass; }
1319 
1320  if((beta==-2) && (zcut==0.1)){
1321  fShapesVar[42]=SymParam;
1322  fShapesVar[43]=DeltaR;
1323  fShapesVar[44]=GroomedPt;
1324  fShapesVar[45]=NGroomedBranches;
1325  fShapesVar[46]=GroomedMass; }
1326 
1327  return;
1328 
1329 
1330 }
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
Bool_t fCentSelectOn
Random number generator.
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:114
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:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
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:112
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:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Double_t Px() const
Definition: AliEmcalJet.h:99
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
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)
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta, Int_t ReclusterAlgo)
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()
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:102
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:44
const AliEmcalPythiaInfo * GetPythiaInfo() const
unsigned short UShort_t
Definition: External.C:28
Double_t Pz() const
Definition: AliEmcalJet.h:101
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
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:256
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t M() const
Definition: AliEmcalJet.h:113
Container for jet within the EMCAL jet framework.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const