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