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