AliPhysics  a8fcd8c (a8fcd8c)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalMissingEnergy.cxx
Go to the documentation of this file.
1 //
2 // N-subjettiness, subjet analysis
3 //
4 // Summer Student: C.Durante - Supervisor: 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 "TMatrixD.h"
20 #include "TMatrixDSym.h"
21 #include "TMatrixDSymEigen.h"
22 #include "TVector3.h"
23 #include "TVector2.h"
24 #include "AliVCluster.h"
25 #include "AliVTrack.h"
26 #include "AliEmcalJet.h"
27 #include "AliRhoParameter.h"
28 #include "AliLog.h"
29 #include "AliEmcalParticle.h"
30 #include "AliMCEvent.h"
31 #include "AliGenPythiaEventHeader.h"
32 #include "AliAODMCHeader.h"
33 #include "AliMCEvent.h"
34 #include "AliAnalysisManager.h"
35 #include "AliJetContainer.h"
36 #include "AliParticleContainer.h"
37 #include "AliEmcalPythiaInfo.h"
38 #include "TRandom3.h"
39 #include "AliPicoTrack.h"
40 #include "AliEmcalJetFinder.h"
41 #include "AliAODEvent.h"
42 #include "AliAnalysisTaskEmcalMissingEnergy.h"
43 #include <algorithm>
44 
45 using std::cout;
46 using std::endl;
47 
48 ClassImp(AliAnalysisTaskEmcalMissingEnergy)
49 
50 //________________________________________________________________________
51 AliAnalysisTaskEmcalMissingEnergy::AliAnalysisTaskEmcalMissingEnergy() :
52 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalMissingEnergy", kTRUE),
53  fContainer(0),
54  fMinFractionShared(0),
55  fJetShapeType(kData),
56  fJetShapeSub(kNoSub),
57  fJetSelection(kInclusive),
58  fSubstructureVar(0),
59  fPtThreshold(-9999.),
60  fRMatching(0.2),
61  fJetRadius(0.4),
62  fSubJetRadius(0.2),
63  fminpTTrig(20.),
64  fmaxpTTrig(50.),
65  fangWindowRecoil(0.6),
66  fSemigoodCorrect(0),
67  fHolePos(0),
68  fHoleWidth(0),
69  fCentSelectOn(kTRUE),
70  fHTon(kTRUE),
71  fCentMin(0),
72  fCentMax(10),
73  fh2ResponseUW(0x0),
74  fh2ResponseW(0x0),
75  fPhiJetCorr6(0x0),
76  fPhiJetCorr7(0x0),
77  fEtaJetCorr6(0x0),
78  fEtaJetCorr7(0x0),
79  fPtJetCorr(0x0),
80  fPtJet(0x0),
81  fhpTjetpT(0x0),
82  fhPt(0x0),
83  fhCentrality(0x0),
84  fhJetPt(0x0),
85  fhJetPhi(0x0),
86  fhTrackPt(0x0),
87  fhTriggerPt(0x0),
88  fhJetTriggeredPt(0x0),
89  fhTriggerVSjetPt(0x0),
90  fhdPhiTrigVSjetPt(0x0),
91  fhdPhiTrigVSjetPtNOCUT(0x0),
92  fhdPhiPartVSjetPt(0x0),
93  fhdPhiPartVSjetPtVSPartPt(0x0),
94  fhPhi(0x0),
95  fSubstructure(0),
96  fHadronTrigger(0),
97  fhTau1(0x0),
98  fhTau2(0x0),
99  fhTau3(0x0),
100  fhTau1vsTau2(0x0),
101  fhTau1vsTau3(0x0),
102  fhTau2vsTau3(0x0),
103  fhTau2OverTau1(0x0),
104  fhTau3OverTau2(0x0),
105  fhTau2vsJetPt(0x0),
106  fhTau2OverTau1vsJetPt(0x0),
107  fhTau1vsTau2vsTau3(0x0),
108  fhTrackPt_JT(0x0),
109  fJetPtMin(10),
110  fTrackPtMin(0.15)
111 {
112  SetMakeGeneralHistograms(kTRUE);
113 }
114 
115 //________________________________________________________________________
116 AliAnalysisTaskEmcalMissingEnergy::AliAnalysisTaskEmcalMissingEnergy(const char *name) :
117  AliAnalysisTaskEmcalJet(name, kTRUE),
118  fContainer(0),
119  fMinFractionShared(0),
120  fJetShapeType(kData),
121  fJetShapeSub(kNoSub),
122  fJetSelection(kInclusive),
123  fSubstructureVar(0),
124  fPtThreshold(-9999.),
125  fRMatching(0.2),
126  fJetRadius(0.4),
127  fSubJetRadius(0.2),
128  fminpTTrig(20.),
129  fmaxpTTrig(50.),
130  fangWindowRecoil(0.6),
131  fSemigoodCorrect(0),
132  fHolePos(0),
133  fHoleWidth(0),
134  fCentSelectOn(kTRUE),
135  fHTon(kTRUE),
136  fCentMin(0),
137  fCentMax(10),
138  fh2ResponseUW(0x0),
139  fh2ResponseW(0x0),
140  fPhiJetCorr6(0x0),
141  fPhiJetCorr7(0x0),
142  fEtaJetCorr6(0x0),
143  fEtaJetCorr7(0x0),
144  fPtJetCorr(0x0),
145  fPtJet(0x0),
146  fhpTjetpT(0x0),
147  fhPt(0x0),
148  fhCentrality(0x0),
149  fhJetPt(0x0),
150  fhJetPhi(0x0),
151  fhTrackPt(0x0),
152  fhTriggerPt(0x0),
153  fhJetTriggeredPt(0x0),
154  fhTriggerVSjetPt(0x0),
155  fhdPhiTrigVSjetPt(0x0),
156  fhdPhiTrigVSjetPtNOCUT(0x0),
157  fhdPhiPartVSjetPt(0x0),
158  fhdPhiPartVSjetPtVSPartPt(0x0),
159  fhPhi(0x0),
160  fSubstructure(0),
161  fHadronTrigger(0),
162  fhTau1(0x0),
163  fhTau2(0x0),
164  fhTau3(0x0),
165  fhTau1vsTau2(0x0),
166  fhTau1vsTau3(0x0),
167  fhTau2vsTau3(0x0),
168  fhTau2OverTau1(0x0),
169  fhTau3OverTau2(0x0),
170  fhTau2vsJetPt(0x0),
171  fhTau2OverTau1vsJetPt(0x0),
172  fhTau1vsTau2vsTau3(0x0),
173  fhTrackPt_JT(0x0),
174  fJetPtMin(10),
175  fTrackPtMin(0.15)
176 {
177  // Standard constructor.
178 
179  SetMakeGeneralHistograms(kTRUE);
180 
181  DefineOutput(1, TTree::Class());
182 
183 }
184 
185 //________________________________________________________________________
186 AliAnalysisTaskEmcalMissingEnergy::~AliAnalysisTaskEmcalMissingEnergy()
187 {
188  // Destructor.
189 }
190 
191 //________________________________________________________________________
192 void AliAnalysisTaskEmcalMissingEnergy::UserCreateOutputObjects()
193 {
194  // Create user output.
195 
196  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
197 
198  // trees
199 
200  fSubstructure = new TTree("fSubstructure","fSubstructure");
201  fHadronTrigger = new TTree("fHadronTrigger","fHadronTrigger");
202 
203  const Int_t nVarS = 6 * 3 + 1;
204  const Int_t nVarH = 5;
205 
206  fSubstructureVar = new Float_t [nVarS];
207  fHadronTriggerVar = new Float_t [nVarH];
208 
209  TString *fSubstructureVarNames = new TString [nVarS];
210  TString *fHadronTriggerVarNames = new TString [nVarH];
211 
212  fSubstructureVarNames[0] = "JetPt";
213  fSubstructureVarNames[1] = "Tau1";
214  fSubstructureVarNames[2] = "Tau2";
215  fSubstructureVarNames[3] = "Tau3";
216  fSubstructureVarNames[4] = "DeltaR2hardest";
217  fSubstructureVarNames[5] = "DeltaPhi2hardest";
218  fSubstructureVarNames[6] = "JetPtMatch";
219  fSubstructureVarNames[7] = "Tau1Match";
220  fSubstructureVarNames[8] = "Tau2Match";
221  fSubstructureVarNames[9] = "Tau3Match";
222  fSubstructureVarNames[10] = "DeltaR2hardestMatch";
223  fSubstructureVarNames[11] = "DeltaPhi2hardestMatch";
224  fSubstructureVarNames[12] = "JetPtDet";
225  fSubstructureVarNames[13] = "Tau1Det";
226  fSubstructureVarNames[14] = "Tau2Det";
227  fSubstructureVarNames[15] = "Tau3Det";
228  fSubstructureVarNames[16] = "DeltaR2hardestDet";
229  fSubstructureVarNames[17] = "DeltaPhi2hardestDet";
230  fSubstructureVarNames[18] = "weightPythia";
231 
232  fHadronTriggerVarNames[0] = "TriggerPt";
233  fHadronTriggerVarNames[1] = "JetPt";
234  fHadronTriggerVarNames[2] = "DeltaPhiJetTrigger";
235  fHadronTriggerVarNames[3] = "PartPt";
236  fHadronTriggerVarNames[4] = "DeltaPhiJetParticles";
237 
238  for(Int_t ivar = 0; ivar<nVarS; ivar++) {
239  cout<<"looping over variables to create the branches of my Substructure Tree"<<endl;
240  fSubstructure->Branch(fSubstructureVarNames[ivar].Data(), &fSubstructureVar[ivar], Form("%s/F", fSubstructureVarNames[ivar].Data()));
241  }
242 
243  for(Int_t ivar = 0; ivar<nVarH; ivar++) {
244  cout<<"looping over variables to create the branches of my HadronTrigger Tree"<<endl;
245  fHadronTrigger->Branch(fHadronTriggerVarNames[ivar].Data(), &fHadronTriggerVar[ivar], Form("%s/F", fHadronTriggerVarNames[ivar].Data()));
246  }
247 
248  // histograms
249 
250  fhCentrality= new TH1F("fhCentrality", "Centrality;c", 200, 0, 100);
251  fOutput->Add(fhCentrality);
252 
253  fhJetPt= new TH1F("fhJetPt", "fhJetPt;Jet p_{T} (GeV)", 100, 0, 100);
254  fOutput->Add(fhJetPt);
255 
256  fhJetPhi= new TH1F("fhJetPhi", "fhJetPhi;Jet #phi", 100, 0, 2*TMath::Pi());
257  fOutput->Add(fhJetPhi);
258 
259  fhTrackPt= new TH1F("fhTrackPt", "fhTrackPt;Track p_{T} (GeV)", 100, 0, 100);
260  fOutput->Add(fhTrackPt);
261 
262  fhTrackPt_JT= new TH1F("fhTrackPt_JT", "fhTrackPt_JT;Track p_{T} (GeV)", 100, 0, 100);
263  fOutput->Add(fhTrackPt_JT);
264 
265  fhTriggerPt= new TH1F("fhTriggerPt", "fhTriggerPt;Trigger p_{T} (GeV)", 100, 0, 100);
266  fOutput->Add(fhTriggerPt);
267 
268  fhJetTriggeredPt= new TH1F("fhJetTriggeredPt", "fhJetTriggeredPt;Triggered Jet p_{T} (GeV)", 100, 0, 100);
269  fOutput->Add(fhJetTriggeredPt);
270 
271  fhTriggerVSjetPt= new TH2F("fhTriggerVSjetPt", "fhTriggerVSjetPt;Trigger p_{T} (GeV);Jet p_{T} (GeV)", 100, 0, 100, 100, 0, 100);
272  fOutput->Add(fhTriggerVSjetPt);
273 
274  fhdPhiTrigVSjetPtNOCUT= new TH2F("fhdPhiTrigVSjetPtNOCUT", "fhdPhiTrigVSjetPtNOCUT;#Delta#phi (Trig-Jet);Jet p_{T} (GeV)", 100, -0.5*TMath::Pi(), 1.5*TMath::Pi(), 100, 0, 100);
275  fOutput->Add(fhdPhiTrigVSjetPtNOCUT);
276 
277  fhdPhiTrigVSjetPt= new TH2F("fhdPhiTrigVSjetPt", "fhdPhiTrigVSjetPt;#Delta#phi (Trig-Jet);Jet p_{T} (GeV)", 100, TMath::Pi() - 0.6, TMath::Pi(), 100, 0, 100);
278  fOutput->Add(fhdPhiTrigVSjetPt);
279 
280  fhdPhiPartVSjetPt= new TH2F("fhdPhiPartVSjetPt", "fhdPhiPartVSjetPt;#Delta#phi (Part-Jet);Jet p_{T} (GeV)", 100, -0.5*TMath::Pi(), 1.5*TMath::Pi(), 100, 0, 100);
281  fOutput->Add(fhdPhiPartVSjetPt);
282 
283  fhdPhiPartVSjetPtVSPartPt= new TH3F("fhdPhiPartVSjetPtVSPartPt", "fhdPhiPartVSjetPtVSPartPt;#Delta#phi (Part-Jet);Jet p_{T} (GeV);Part p_{T} (GeV)", 100, -0.5*TMath::Pi(), 1.5*TMath::Pi(), 100, 0, 100, 100, 0, 100);
284  fOutput->Add(fhdPhiPartVSjetPtVSPartPt);
285 
286  // substructure output
287 
288  fhTau1 = new TH1F("fhTau1", "fhTau1;#tau_{1}", 100, 0, 1);
289  fOutput->Add(fhTau1);
290 
291  fhTau2 = new TH1F("fhTau2", "fhTau2;#tau_{2}", 100, 0, 1);
292  fOutput->Add(fhTau2);
293 
294  fhTau3 = new TH1F("fhTau3", "fhTau3;#tau_{3}", 100, 0, 1);
295  fOutput->Add(fhTau3);
296 
297  fhTau1vsTau2 = new TH2F("fhTau1vsTau2", "fhTau1vsTau2;#tau_{1};#tau_{2}", 100, 0, 1, 100, 0, 1);
298  fOutput->Add(fhTau1vsTau2);
299 
300  fhTau1vsTau3 = new TH2F("fhTau1vsTau3", "fhTau1vsTau3;#tau_{1};#tau_{3}", 100, 0, 1, 100, 0, 1);
301  fOutput->Add(fhTau1vsTau3);
302 
303  fhTau2vsTau3 = new TH2F("fhTau2vsTau3", "fhTau2vsTau3;#tau_{2};#tau_{3}", 100, 0, 1, 100, 0, 1);
304  fOutput->Add(fhTau2vsTau3);
305 
306  fhTau1vsTau2vsTau3 = new TH3F("fhTau1vsTau2vsTau3", "fhTau1vsTau2vsTau3;#tau_{1};#tau_{2};#tau_{3}", 100, 0, 1, 100, 0, 1, 100, 0, 1);
307  fOutput->Add(fhTau1vsTau2vsTau3);
308 
309  fhTau2OverTau1 = new TH1F("fhTau2OverTau1", "fhTau2OverTau1;#tau_{2}/#tau_{1}", 100, 0, 1);
310  fOutput->Add(fhTau2OverTau1);
311 
312  fhTau3OverTau2 = new TH1F("fhTau3OverTau2", "fhTau3OverTau2;#tau_{3}/#tau_{2}", 100, 0, 1);
313  fOutput->Add(fhTau3OverTau2);
314 
315 
316  fhTau2vsJetPt = new TH2F("fhTau2vsJetPt", "fhTau2vsJetPt;#tau_{2}; Jet p_{T} (GeV)", 100, 0, 1, 100, 0, 100);
317  fOutput->Add(fhTau2vsJetPt);
318 
319  fhTau2OverTau1vsJetPt = new TH2F("fhTau2OverTau1vsJetPt", "fhTau2OverTau1vsJetPt;#tau_{2}/#tau_{1}; Jet p_{T} (GeV)", 100, 0, 1, 100, 0, 100);
320  fOutput->Add(fhTau2OverTau1vsJetPt);
321 
322  fOutput->Add(fSubstructure);
323  fOutput->Add(fHadronTrigger);
324  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
325 
326 }
327 
328 //________________________________________________________________________
329 Bool_t AliAnalysisTaskEmcalMissingEnergy::Run()
330 {
331  // Run analysis code here, if needed. It will be executed before FillHistograms().
332 
333  return kTRUE;
334 }
335 
336 //________________________________________________________________________
337 Bool_t AliAnalysisTaskEmcalMissingEnergy::FillHistograms()
338 {
339 
340  if (fCentSelectOn) {
341  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
342  }
343  fhCentrality->Fill(fCent);
344 
345  const Double_t eta_max = 0.9;
346  const Double_t phi_h = 0.6;
347 
348  // const Double_t epsilon = 0.00001;
349 
350  AliParticleContainer *partCont = GetParticleContainer(0);
351  TClonesArray *tracksArray = partCont->GetArray();
352 
353  Float_t kWeight=1;
354 
355  if(!partCont || !tracksArray) return 0;
356  AliPicoTrack *trigger = 0x0;
357  AliPicoTrack *track = 0x0;
358 
359  AliEmcalJet* jetTEMP = NULL;
360  AliJetContainer *JetCont = GetJetContainer(0);
361 
362  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){ //filling track histogram
363  track = static_cast<AliPicoTrack*>(tracksArray->At(iTrack));
364  if (!track || track->Pt()<1) continue;
365  fhTrackPt->Fill(track->Pt());
366  }
367 
368  if(JetCont) {
369  JetCont->ResetCurrentID();
370  while((jetTEMP = JetCont->GetNextAcceptJet())) {
371  if (!jetTEMP) continue;
372  if (jetTEMP->Pt() < fJetPtMin) continue;
373  if (TMath::Abs(jetTEMP->Eta()) > eta_max - fJetRadius) continue;
374  fhJetPt->Fill(jetTEMP->Pt());
375  fhJetPhi->Fill(jetTEMP->Phi());
376  }
377  }
378 
379  // start Hadron Trigger
380 
381  Int_t trigger_index = SelectTrigger(fminpTTrig,fmaxpTTrig);
382 
383  if(trigger_index >= 0 && fHTon) {
384  trigger = static_cast<AliPicoTrack*>(tracksArray->At(trigger_index));
385  if (trigger) {
386  fhTriggerPt->Fill(trigger->Pt());
387 
388  if(JetCont) {
389  JetCont->ResetCurrentID();
390  while((jetTEMP = JetCont->GetNextAcceptJet())) {
391  if (!jetTEMP) continue;
392  if (jetTEMP->Pt() < fJetPtMin) continue;
393  if (TMath::Abs(jetTEMP->Eta()) > eta_max - fJetRadius) continue;
394 
395  Float_t dphi = RelativePhi(trigger->Phi(), jetTEMP->Phi()); // -pi to pi
396  Float_t dphiFancy = RelativePhiFancy(trigger->Phi(), jetTEMP->Phi()); // -0.5pi to 1.5pi
397 
398  fhdPhiTrigVSjetPtNOCUT->Fill(dphiFancy,jetTEMP->Pt());
399 
400  if (TMath::Abs(dphi) < TMath::Pi() - phi_h) continue; // trigger selection
401 
402  fhJetTriggeredPt->Fill(jetTEMP->Pt());
403  fhTriggerVSjetPt->Fill(trigger->Pt(),jetTEMP->Pt());
404  fhdPhiTrigVSjetPt->Fill(TMath::Abs(dphi),jetTEMP->Pt());
405 
406  //now I have a trigger and a recoil jet. Let's loop over all particles!
407  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
408  track = static_cast<AliPicoTrack*>(tracksArray->At(iTrack));
409  if (!track) continue;
410 
411  if(TMath::Abs(track->Eta())>eta_max) continue;
412  if(track->Pt()<fTrackPtMin) continue;
413  //if(track->GetTrackType() == 2) continue;
414 
415  Float_t DPHI = RelativePhiFancy(jetTEMP->Phi(), track->Phi()); // -0.5 pi to 1.5 pi
416  fhdPhiPartVSjetPt->Fill(DPHI,jetTEMP->Pt());
417  fhdPhiPartVSjetPtVSPartPt->Fill(DPHI,jetTEMP->Pt(),track->Pt());
418  fhTrackPt_JT->Fill(track->Pt());
419 
420  fHadronTriggerVar[0] = trigger->Pt();
421  fHadronTriggerVar[1] = jetTEMP->Pt();
422  fHadronTriggerVar[2] = TMath::Abs(dphi);
423  fHadronTriggerVar[3] = track->Pt();
424  fHadronTriggerVar[4] = DPHI;
425 
426  fHadronTrigger->Fill();
427  }
428  }
429  }
430  }
431  }
433 
434 
436 
437  Float_t *shape, *shapePart, *shapeDet;
438  AliEmcalJet *jet1 = NULL;
439 
440  if(JetCont){
441  JetCont->ResetCurrentID();
442 
443  while((jet1 = JetCont->GetNextAcceptJet())) { // loop over main Jets found in the event
444  if(!jet1) continue;
445  if (jet1->Pt() < fJetPtMin) continue;
446  if (TMath::Abs(jet1->Eta()) > eta_max - fJetRadius) continue;
447 
448  AliEmcalJet* jet2 = 0x0;
449  AliEmcalJet* jet3 = 0x0;
450  AliEmcalJet* jetUS = 0x0;
451  Int_t ifound=0;
452  Int_t ilab=-1;
453  if (fJetShapeType != kData) {
454  const AliEmcalPythiaInfo *partonsInfo = 0x0;
455  AliJetContainer *jetContTrue = GetJetContainer(1);
456  AliJetContainer *jetContUS = GetJetContainer(2);
457 
458  if (fJetShapeSub == kConstSub) {
459  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
460  jetUS = jetContUS->GetJet(i);
461  if(jetUS->GetLabel()==jet1->GetLabel()) {
462  ifound++;
463  if(ifound==1) ilab = i;
464  }
465  }
466  if(ilab == -1) continue;
467  jetUS=jetContUS->GetJet(ilab);
468  jet2=jetUS->ClosestJet();
469  }
470 
471  if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet();
472  if (!jet2) {
473  Printf("jet2 does not exist, returning");
474  continue;
475  }
476 
477  if(fJetShapeType==kDetEmbPart){
478  AliJetContainer *jetContPart=GetJetContainer(3);
479  jet3=jet2->ClosestJet();
480 
481  if(!jet3){
482  Printf("jet3 does not exist, returning");
483  continue;
484  }
485  cout<<"jet 3 exists"<<jet3->Pt()<<endl;
486 
487  Double_t fraction=0;
488  if(!(fJetShapeSub==kConstSub)) fraction = JetCont->GetFractionSharedPt(jet1);
489  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
490  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
491 
492  if(fraction<fMinFractionShared) continue;
493 
494  partonsInfo = GetPythiaInfo();
495  if(!partonsInfo) return 0;
496  kWeight = partonsInfo->GetPythiaEventWeight();
497 
498  }
499 
500  }
501 
502  shape = N_subjettiness(jet1,0);
503  if (shape[6]<5.) continue; // check that the tau computation was correct, no errors
504  for (Int_t i=0; i<6; i++) {
505  fSubstructureVar[i] = shape[i];
506  }
507 
508  if (fJetShapeType == kData) {
509  for (Int_t i=6; i<18; i++) {
510  fSubstructureVar[i] = -2;
511  }
512 
513  fSubstructureVar[18]= 1;
514  }
515  else if(fJetShapeType==kDetEmbPart){
516  shapePart = N_subjettiness(jet3,3); // particle level
517  shapeDet = N_subjettiness(jet2,2); // detector level
518  if (shapePart[6]<5. || shapeDet[6]<5.) continue; // check that the tau computation was correct, no errors
519  for (Int_t i=0; i<6; i++) {
520  fSubstructureVar[i+6] = shapePart[i];
521  fSubstructureVar[i+12] = shapeDet[i];
522  }
523 
524  fSubstructureVar[18]= kWeight;
525  }
526 
527  if (shape[1]>=0. && shape[2]<0. && shape[3]<0.) {
528  fhTau1->Fill(shape[1]);
529  }
530  else if (shape[1]>=0. && shape[2]>0. && shape[3]<0.) {
531  fhTau1->Fill(shape[1]);
532  fhTau2->Fill(shape[2]);
533  fhTau1vsTau2->Fill(shape[1],shape[2]);
534  fhTau2OverTau1->Fill(shape[2]/shape[1]);
535  fhTau2vsJetPt->Fill(shape[2],shape[0]);
536  fhTau2OverTau1vsJetPt->Fill(shape[2]/shape[1],shape[0]);
537  }
538  else if (shape[1]>=0. && shape[2]>0. && shape[3]>0.) {
539  fhTau1->Fill(shape[1]);
540  fhTau2->Fill(shape[2]);
541  fhTau3->Fill(shape[3]);
542  fhTau1vsTau2->Fill(shape[1],shape[2]);
543  fhTau1vsTau3->Fill(shape[1],shape[3]);
544  fhTau2vsTau3->Fill(shape[2],shape[3]);
545  fhTau1vsTau2vsTau3->Fill(shape[1],shape[2],shape[3]);
546  fhTau2OverTau1->Fill(shape[2]/shape[1]);
547  fhTau3OverTau2->Fill(shape[3]/shape[2]);
548  fhTau2vsJetPt->Fill(shape[2],shape[0]);
549  fhTau2OverTau1vsJetPt->Fill(shape[2]/shape[1],shape[0]);
550  }
551 
552  fSubstructure->Fill();
553  }
554  }
555 
556 
558 
559  return kTRUE;
560 }
561 
562 //________________________________________________________________________
563 Float_t AliAnalysisTaskEmcalMissingEnergy::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0) {
564  //calc subtracted jet mass
565  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
566  return jet->GetShapeProperties()->GetSecondOrderSubtracted();
567  else
568  return jet->M();
569 }
570 
571 //________________________________________________________________________
572 Float_t AliAnalysisTaskEmcalMissingEnergy::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){
573 
574  AliJetContainer *jetCont = GetJetContainer(jetContNb);
575  if (!jet->GetNumberOfTracks())
576  return 0;
577  Double_t den=0.;
578  Double_t num = 0.;
579  AliVParticle *vp1 = 0x0;
580  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
581  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
582 
583  if (!vp1){
584  Printf("AliVParticle associated to constituent not found");
585  continue;
586  }
587 
588  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
589  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
590  Double_t dr = TMath::Sqrt(dr2);
591  num=num+vp1->Pt()*dr;
592  den=den+vp1->Pt();
593  }
594  return num/den;
595 }
596 
597 //________________________________________________________________________
598 Float_t AliAnalysisTaskEmcalMissingEnergy::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) {
599 
600  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
601  return jet->GetShapeProperties()->GetSecondOrderSubtractedAngularity();
602  else
603  return Angularity(jet, jetContNb);
604 
605 }
606 
607 
608 //________________________________________________________________________
609 Float_t AliAnalysisTaskEmcalMissingEnergy::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){
610 
611  AliJetContainer *jetCont = GetJetContainer(jetContNb);
612  if (!jet->GetNumberOfTracks())
613  return 0;
614  Double_t den = 0.;
615  Double_t num = 0.;
616  AliVParticle *vp1 = 0x0;
617  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
618  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
619 
620  if (!vp1){
621  Printf("AliVParticle associated to constituent not found");
622  continue;
623  }
624 
625  num=num+vp1->Pt()*vp1->Pt();
626  den=den+vp1->Pt();
627  }
628  return TMath::Sqrt(num)/den;
629 }
630 
631 //________________________________________________________________________
632 Float_t AliAnalysisTaskEmcalMissingEnergy::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) {
633  //calc subtracted jet mass
634  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
635  return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
636  else
637  return PTD(jet, jetContNb);
638 
639 }
640 
641 //_____________________________________________________________________________
642 Float_t AliAnalysisTaskEmcalMissingEnergy::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){
643 
644  AliJetContainer *jetCont = GetJetContainer(jetContNb);
645  if (!jet->GetNumberOfTracks())
646  return 0;
647  Double_t mxx = 0.;
648  Double_t myy = 0.;
649  Double_t mxy = 0.;
650  int nc = 0;
651  Double_t sump2 = 0.;
652  Double_t pxjet=jet->Px();
653  Double_t pyjet=jet->Py();
654  Double_t pzjet=jet->Pz();
655 
656 
657  //2 general normalized vectors perpendicular to the jet
658  TVector3 ppJ1(pxjet, pyjet, pzjet);
659  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
660  ppJ3.SetMag(1.);
661  TVector3 ppJ2(-pyjet, pxjet, 0);
662  ppJ2.SetMag(1.);
663  AliVParticle *vp1 = 0x0;
664  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
665  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
666 
667  if (!vp1){
668  Printf("AliVParticle associated to constituent not found");
669  continue;
670  }
671 
672  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
673 
674  //local frame
675  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
676  TVector3 pPerp = pp - pLong;
677  //projection onto the two perpendicular vectors defined above
678 
679  Float_t ppjX = pPerp.Dot(ppJ2);
680  Float_t ppjY = pPerp.Dot(ppJ3);
681  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
682  if(ppjT<=0) return 0;
683 
684  mxx += (ppjX * ppjX / ppjT);
685  myy += (ppjY * ppjY / ppjT);
686  mxy += (ppjX * ppjY / ppjT);
687  nc++;
688  sump2 += ppjT;}
689 
690  if(nc<2) return 0;
691  if(sump2==0) return 0;
692  // Sphericity Matrix
693  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
694  TMatrixDSym m0(2,ele);
695 
696  // Find eigenvectors
697  TMatrixDSymEigen m(m0);
698  TVectorD eval(2);
699  TMatrixD evecm = m.GetEigenVectors();
700  eval = m.GetEigenValues();
701  // Largest eigenvector
702  int jev = 0;
703  // cout<<eval[0]<<" "<<eval[1]<<endl;
704  if (eval[0] < eval[1]) jev = 1;
705  TVectorD evec0(2);
706  // Principle axis
707  evec0 = TMatrixDColumn(evecm, jev);
708  Double_t compx=evec0[0];
709  Double_t compy=evec0[1];
710  TVector2 evec(compx, compy);
711  Double_t circ=0;
712  if(jev==1) circ=2*eval[0];
713  if(jev==0) circ=2*eval[1];
714 
715  return circ;
716 }
717 
718 
719 //________________________________________________________________________
720 Float_t AliAnalysisTaskEmcalMissingEnergy::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) {
721  //calc subtracted jet mass
722 
723  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
724  return jet->GetShapeProperties()->GetSecondOrderSubtractedCircularity();
725  else
726  return Circularity(jet, jetContNb);
727 }
728 
729 //________________________________________________________________________
730 Float_t AliAnalysisTaskEmcalMissingEnergy::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){
731 
732  AliJetContainer *jetCont = GetJetContainer(jetContNb);
733  if (!jet->GetNumberOfTracks())
734  return 0;
735  Double_t den=0.;
736  Double_t num = 0.;
737  AliVParticle *vp1 = 0x0;
738  AliVParticle *vp2 = 0x0;
739  std::vector<int> ordindex;
740  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
741  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
742  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
743 
744  if(ordindex.size()<2) return -1;
745 
746  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
747  if (!vp1){
748  Printf("AliVParticle associated to Leading constituent not found");
749  return -1;
750  }
751 
752  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
753  if (!vp2){
754  Printf("AliVParticle associated to Subleading constituent not found");
755  return -1;
756  }
757 
758 
759  num=vp1->Pt();
760  den=vp2->Pt();
761  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
762 
763  return num-den;
764 }
765 
766 //________________________________________________________________________
767 Float_t AliAnalysisTaskEmcalMissingEnergy::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) {
768  //calc subtracted jet mass
769 
770  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
771  return jet->GetShapeProperties()->GetSecondOrderSubtractedLeSub();
772  else
773  return LeSub(jet, jetContNb);
774 
775 }
776 
777 //________________________________________________________________________
778 Float_t AliAnalysisTaskEmcalMissingEnergy::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0) {
779  //calc subtracted jet mass
780 
781  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
782  return jet->GetShapeProperties()->GetSecondOrderSubtractedConstituent();
783  else
784  return jet->GetNumberOfTracks();
785 
786 }
787 
788 
789 //______________________________________________________________________________
790 Float_t AliAnalysisTaskEmcalMissingEnergy::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){
791 
792  AliJetContainer *jetCont = GetJetContainer(jetContNb);
793  if (!jet->GetNumberOfTracks())
794  return 0;
795  Double_t mxx = 0.;
796  Double_t myy = 0.;
797  Double_t mxy = 0.;
798  int nc = 0;
799  Double_t sump2 = 0.;
800 
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 ppt=vp1->Pt();
811  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
812 
813  Double_t deta = vp1->Eta()-jet->Eta();
814  mxx += ppt*ppt*deta*deta;
815  myy += ppt*ppt*dphi*dphi;
816  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
817  nc++;
818  sump2 += ppt*ppt;
819 
820  }
821  if(nc<2) return 0;
822  if(sump2==0) return 0;
823  // Sphericity Matrix
824  Double_t ele[4] = {mxx , mxy , mxy , myy };
825  TMatrixDSym m0(2,ele);
826 
827  // Find eigenvectors
828  TMatrixDSymEigen m(m0);
829  TVectorD eval(2);
830  TMatrixD evecm = m.GetEigenVectors();
831  eval = m.GetEigenValues();
832  // Largest eigenvector
833  int jev = 0;
834  // cout<<eval[0]<<" "<<eval[1]<<endl;
835  if (eval[0] < eval[1]) jev = 1;
836  TVectorD evec0(2);
837  // Principle axis
838  evec0 = TMatrixDColumn(evecm, jev);
839  Double_t compx=evec0[0];
840  Double_t compy=evec0[1];
841  TVector2 evec(compx, compy);
842  Double_t sig=0;
843  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
844  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
845 
846  return sig;
847 
848 }
849 
850 //________________________________________________________________________
851 Float_t AliAnalysisTaskEmcalMissingEnergy::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) {
852  //calc subtracted jet mass
853 
854  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
855  return jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2();
856  else
857  return Sigma2(jet, jetContNb);
858 
859 }
860 
861 //________________________________________________________________________
862 Int_t AliAnalysisTaskEmcalMissingEnergy::SelectTrigger(Float_t minpT, Float_t maxpT){
863 
864  AliParticleContainer *partCont = GetParticleContainer(0);
865  TClonesArray *tracksArray = partCont->GetArray();
866 
867  if(!partCont || !tracksArray) return -99999;
868  AliAODTrack *track = 0x0;
869  AliEmcalParticle *emcPart = 0x0;
870  AliPicoTrack *picoTrack = 0x0;
871 
872  TList *trackList = new TList();
873  Int_t triggers[100];
874  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
875  Int_t iTT = 0;
876 
877  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
878 
879  if ((fJetShapeSub == kNoSub) || (fJetShapeSub == kDerivSub)) {
880  picoTrack = static_cast<AliPicoTrack*>(tracksArray->At(iTrack));
881  if (!picoTrack) continue;
882 
883 
884  if(TMath::Abs(picoTrack->Eta())>0.9) continue;
885  if(picoTrack->Pt()<0.15) continue;
886  if(picoTrack->GetTrackType() == 2) continue;
887 
888  //if ((picoTrack->Pt()>8) && (picoTrack->Pt()<9)) Printf("picoTrackLabel = %d", picoTrack->GetTrackType());
889 
890  if ((picoTrack->Pt() >= minpT) && (picoTrack->Pt()< maxpT)) {
891  trackList->Add(picoTrack);
892  triggers[iTT] = iTrack;
893  iTT++;
894  }
895  }
896  else if (fJetShapeSub == kConstSub){
897  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
898  if (!emcPart) continue;
899  if(TMath::Abs(emcPart->Eta())>0.9) continue;
900  if (emcPart->Pt()<0.15) continue;
901 
902  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
903  trackList->Add(emcPart);
904  triggers[iTT] = iTrack;
905  iTT++;
906  }
907  }
908  else{
909  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
910  if (!track) continue;
911  if(TMath::Abs(track->Eta())>0.9) continue;
912  if (track->Pt()<0.15) continue;
913  if (!(track->TestFilterBit(768))) continue;
914 
915  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
916  trackList->Add(track);
917  triggers[iTT] = iTrack;
918  iTT++;
919 
920  }
921  }
922  }
923 
924  if (iTT == 0) return -99999;
925  Int_t nbRn = 0, index = 0 ;
926  TRandom3* random = new TRandom3(0);
927  nbRn = random->Integer(iTT);
928 
929  index = triggers[nbRn];
930  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
931  return index;
932 
933 }
934 
935 //__________________________________________________________________________________
936 Double_t AliAnalysisTaskEmcalMissingEnergy::RelativePhi(Double_t mphi,Double_t vphi){
937 
938  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
939  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
940  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
941  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
942  double dphi = mphi-vphi;
943  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
944  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
945  return dphi;//dphi in [-Pi, Pi]
946 }
947 
948 
949 //________________________________________________________________________
950 Bool_t AliAnalysisTaskEmcalMissingEnergy::RetrieveEventObjects() {
951  //
952  // retrieve event objects
953  //
954  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
955  return kFALSE;
956 
957  return kTRUE;
958 }
959 
960 //_______________________________________________________________________
961 void AliAnalysisTaskEmcalMissingEnergy::Terminate(Option_t *)
962 {
963  // Called once at the end of the analysis.
964 
965  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
966  // if (!fTreeObservableTagging){
967  // Printf("ERROR: fTreeObservable not available");
968  // return;
969  // }
970 
971 }
972 
973 //_______________________________________________________________________
974 Double_t AliAnalysisTaskEmcalMissingEnergy::Minimum(Double_t x, Double_t y)
975 {
976  if (x <= y) return x;
977  else return y;
978 }
979 
980 //_______________________________________________________________________
981 Double_t AliAnalysisTaskEmcalMissingEnergy::Minimum(Double_t x, Double_t y, Double_t z)
982 {
983  Double_t minimum = x;
984  if (y < minimum) minimum = y;
985  if (z < minimum) minimum = z;
986  return minimum;
987 }
988 
989 //_______________________________________________________________________
990 Double_t AliAnalysisTaskEmcalMissingEnergy::R_distance(AliEmcalJet *jet1, AliEmcalJet *jet2)
991 {
992  Double_t dPhi, dEta;
993  dPhi = RelativePhi(jet1->Phi(),jet2->Phi());
994  dEta = jet1->Eta() - jet2->Eta();
995  return TMath::Sqrt(dEta*dEta + dPhi*dPhi);
996 }
997 
998 //_______________________________________________________________________
999 Double_t AliAnalysisTaskEmcalMissingEnergy::R_distance(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
1000 {
1001  Double_t dPhi, dEta;
1002  dPhi = RelativePhi(phi1,phi2);
1003  dEta = eta1 - eta2;
1004  return TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1005 }
1006 
1007 //_______________________________________________________________________
1008 Double_t AliAnalysisTaskEmcalMissingEnergy::R_distance(AliEmcalJet *jet, AliPicoTrack *part)
1009 {
1010  Double_t dPhi, dEta;
1011  dPhi = RelativePhi(jet->Phi(),part->Phi());
1012  dEta = jet->Eta() - part->Eta();
1013  return TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1014 }
1015 
1016 //__________________________________________________________________________________
1017 Double_t AliAnalysisTaskEmcalMissingEnergy::RelativePhiFancy(Double_t mphi,Double_t vphi){
1018 
1019  Double_t dPhi = mphi - vphi;
1020 
1021  if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
1022  if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
1023  if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1024  if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1025 
1026  return dPhi; // [-0.5 pi, 1.5 pi]
1027 }
1028 
1029 
1030 //__________________________________________________________________________________
1031 Int_t * AliAnalysisTaskEmcalMissingEnergy::JetHard(AliEmcalJetFinder *finder){
1032 
1033  Int_t *indeces;
1034  Int_t NumberOfJets = finder->GetNumberOfJets();
1035  AliEmcalJet *temp;
1036  Float_t mom;
1037  Int_t ind;
1038 
1039  Float_t *PtArray;
1040  PtArray = new Float_t[NumberOfJets];
1041  indeces = new Int_t[NumberOfJets];
1042 
1043  for(Int_t i=0;i<NumberOfJets;i++){
1044  temp = finder->GetJet(i);
1045  PtArray[i] = temp->Pt();
1046  indeces[i] = i;
1047  }
1048 
1049  for(Int_t i=0; i<NumberOfJets - 1; i++){
1050  for(Int_t j=0; j<(NumberOfJets - i - 1); j++){
1051  if(PtArray[j]>PtArray[j+1]){
1052  mom=PtArray[j];
1053  PtArray[j]=PtArray[j+1];
1054  PtArray[j+1]=mom;
1055 
1056  ind = indeces[j];
1057  indeces[j]=indeces[j+1];
1058  indeces[j+1]=ind;
1059  }
1060  }
1061  }
1062 
1063  return indeces;
1064 }
1065 
1066 //__________________________________________________________________________________
1067 Float_t AliAnalysisTaskEmcalMissingEnergy::TauDen(AliEmcalJet *mainJet, Int_t jetContNb = 0){
1068 
1069  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1070  AliPicoTrack *mainJetParticle = 0x0;
1071  if (!mainJet->GetNumberOfTracks()) return 0;
1072  Float_t tau_den = 0.;
1073  for (Int_t i = 0; i < mainJet->GetNumberOfTracks(); i++) { //compute tau denominator once and for all
1074  mainJetParticle = static_cast<AliPicoTrack*>(mainJet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1075  tau_den += mainJetParticle->Pt() * fJetRadius;
1076  }
1077  return tau_den;
1078 }
1079 
1080 //__________________________________________________________________________________
1081 Float_t AliAnalysisTaskEmcalMissingEnergy::Tau1Num(AliEmcalJet *mainJet, AliEmcalJet *subJet1hardest, Int_t jetContNb = 0){
1082 
1083  Float_t dR1;
1084  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1085  AliPicoTrack *mainJetParticle = 0x0;
1086  if (!mainJet->GetNumberOfTracks()) return 0;
1087  Float_t tau1_num = 0.;
1088  for (Int_t i = 0; i < mainJet->GetNumberOfTracks(); i++) { // tau1 computation
1089  mainJetParticle = static_cast<AliPicoTrack*>(mainJet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1090  dR1 = R_distance(subJet1hardest,mainJetParticle);
1091  tau1_num += mainJetParticle->Pt() * dR1;
1092  }
1093  return tau1_num;
1094 }
1095 
1096 //__________________________________________________________________________________
1097 Float_t AliAnalysisTaskEmcalMissingEnergy::Tau1Num_full(AliEmcalJet *mainJet, AliEmcalJetFinder *finder, Int_t jetContNb = 0){
1098 
1099  Int_t NumberOfSubjets = finder->GetNumberOfJets();
1100  Float_t dR1, tau1_num_min = 10;
1101  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1102  AliPicoTrack *mainJetParticle = 0x0;
1103  AliEmcalJet *subJet = NULL;
1104  if (!mainJet->GetNumberOfTracks()) return 0.;
1105 
1106  for (Int_t n=0; n<NumberOfSubjets; n++) {
1107  Float_t tau1_num = 0.;
1108  subJet = finder->GetJet(n);
1109  for (Int_t i = 0; i < mainJet->GetNumberOfTracks(); i++) { // tau1 computation
1110  mainJetParticle = static_cast<AliPicoTrack*>(mainJet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1111  dR1 = R_distance(subJet,mainJetParticle);
1112  tau1_num += mainJetParticle->Pt() * dR1;
1113  }
1114  if (tau1_num < tau1_num_min) tau1_num_min = tau1_num;
1115  }
1116  return tau1_num_min;
1117 }
1118 
1119 
1120 //__________________________________________________________________________________
1121 Float_t AliAnalysisTaskEmcalMissingEnergy::Tau2Num(AliEmcalJet *mainJet, AliEmcalJet *subJet1hardest, AliEmcalJet *subJet2hardest, Int_t jetContNb = 0){
1122 
1123  Float_t dR1, dR2, dRmin;
1124  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1125  AliPicoTrack *mainJetParticle = 0x0;
1126  if (!mainJet->GetNumberOfTracks()) return 0;
1127  Float_t tau2_num = 0.;
1128  for (Int_t i = 0; i < mainJet->GetNumberOfTracks(); i++) { //tau2 computation
1129  mainJetParticle = static_cast<AliPicoTrack*>(mainJet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1130  dR1 = R_distance(subJet1hardest,mainJetParticle);
1131  dR2 = R_distance(subJet2hardest,mainJetParticle);
1132  dRmin = Minimum(dR1, dR2);
1133  tau2_num += mainJetParticle->Pt() * dRmin;
1134  }
1135  return tau2_num;
1136 }
1137 
1138 //__________________________________________________________________________________
1139 Float_t AliAnalysisTaskEmcalMissingEnergy::Tau3Num(AliEmcalJet *mainJet, AliEmcalJet *subJet1hardest, AliEmcalJet *subJet2hardest, AliEmcalJet *subJet3hardest, Int_t jetContNb = 0){
1140 
1141  Float_t dR1, dR2, dR3, dRmin;
1142  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1143  AliPicoTrack *mainJetParticle = 0x0;
1144  if (!mainJet->GetNumberOfTracks()) return 0;
1145  Float_t tau3_num = 0.;
1146  for (Int_t i = 0; i < mainJet->GetNumberOfTracks(); i++) { //tau3 computation
1147  mainJetParticle = static_cast<AliPicoTrack*>(mainJet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1148  dR1 = R_distance(subJet1hardest,mainJetParticle);
1149  dR2 = R_distance(subJet2hardest,mainJetParticle);
1150  dR3 = R_distance(subJet3hardest,mainJetParticle);
1151  dRmin = Minimum(dR1,dR2,dR3);
1152  tau3_num += mainJetParticle->Pt() * dRmin;
1153  }
1154  return tau3_num;
1155 }
1156 
1157 
1158 //__________________________________________________________________________________
1159 Float_t* AliAnalysisTaskEmcalMissingEnergy::N_subjettiness(AliEmcalJet *mainJet, Int_t jetContNb = 0){
1160  Float_t *shape; // jetPt, tau1, tau2, tau3, deltaR2hardest, deltaPhi2hardest, control [if (control<5) continue;]
1161  shape = new Float_t[7];
1162  AliJetContainer * jetCont = GetJetContainer(jetContNb);
1163  if (!mainJet->GetNumberOfTracks()) return 0;
1164  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder("SubJetFinder");
1165  Int_t *ordered;
1166  Int_t SubJetCounter;
1167  Float_t tau1_num = 0., tau2_num = 0., tau3_num = 0., tau_den = 0., dR2hardest = 0., dPhi2hardest = 0.;
1168 
1169  AliEmcalJet *subJet = NULL, *subJet1hardest = NULL, *subJet2hardest = NULL, *subJet3hardest = NULL;
1170 
1171  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1172  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1173 
1174  Reclusterer->SetRadius(fSubJetRadius);
1175  Reclusterer->SetJetMinPt(fTrackPtMin);
1176  Reclusterer->SetJetAlgorithm(0); // 0 for anti-kt 1 for kt
1177 
1178  do {
1179  if(!(Reclusterer->AliEmcalJetFinder::Filter(mainJet, jetCont, dVtx))) continue; //reclustering mainJets using the jetfinderobject Reclusterer
1180  SubJetCounter = Reclusterer->GetNumberOfJets(); // Number of reclustered SubJets in each original jet
1181  if(0 >= SubJetCounter) continue;
1182  ordered = JetHard(Reclusterer); // subjet pt ordering to get the hardests
1183  if(NULL == ordered) continue;
1184  for (Int_t i=0; i<SubJetCounter; i++){ //Loops through each SubJet in a reclustered jet
1185  subJet = Reclusterer->GetJet(ordered[i]);
1186  cout<<"SubJet: "<<i+1<<"/"<<SubJetCounter<<" | subJetPt/JetPt: "<<subJet->Pt()<<"/"<<mainJet->Pt()<<endl;
1187  }
1188  tau_den = TauDen(mainJet, jetContNb);
1189  if (tau_den <= 0.) continue;
1190  shape[0] = mainJet->Pt();
1191  subJet1hardest = Reclusterer->GetJet(ordered[SubJetCounter - 1]);
1192  if (NULL == subJet1hardest) continue;
1193  tau1_num = Tau1Num(mainJet,subJet1hardest,jetContNb);
1194 
1195  if (1 == SubJetCounter) {
1196  if (tau1_num == 0.) continue;
1197  shape[1] = tau1_num/tau_den;
1198  shape[2] = -1;
1199  shape[3] = -1;
1200  shape[4] = -1;
1201  shape[5] = -1;
1202  }
1203  else if (2 == SubJetCounter) {
1204  subJet2hardest = Reclusterer->GetJet(ordered[SubJetCounter - 2]);
1205  if (NULL == subJet2hardest) continue;
1206  tau2_num = Tau2Num(mainJet,subJet1hardest,subJet2hardest,jetContNb);
1207  dR2hardest = R_distance(subJet1hardest,subJet2hardest);
1208  dPhi2hardest = TMath::Abs(RelativePhi(subJet1hardest->Phi(),subJet2hardest->Phi()));
1209  if (tau1_num == 0. || tau2_num == 0.) continue;
1210  shape[1] = tau1_num/tau_den;
1211  shape[2] = tau2_num/tau_den;
1212  shape[3] = -1;
1213  shape[4] = dR2hardest;
1214  shape[5] = dPhi2hardest;
1215  }
1216  else if (SubJetCounter > 2) {
1217  subJet2hardest = Reclusterer->GetJet(ordered[SubJetCounter - 2]);
1218  if(NULL == subJet2hardest) continue;
1219  subJet3hardest = Reclusterer->GetJet(ordered[SubJetCounter - 3]);
1220  if(NULL == subJet3hardest) continue;
1221  tau2_num = Tau2Num(mainJet,subJet1hardest,subJet2hardest,jetContNb);
1222  tau3_num = Tau3Num(mainJet,subJet1hardest,subJet2hardest,subJet3hardest,jetContNb);
1223  dR2hardest = R_distance(subJet1hardest,subJet2hardest);
1224  dPhi2hardest = TMath::Abs(RelativePhi(subJet1hardest->Phi(),subJet2hardest->Phi()));
1225  if (tau1_num == 0. || tau2_num == 0. || tau3_num == 0.) continue;
1226  shape[1] = tau1_num/tau_den;
1227  shape[2] = tau2_num/tau_den;
1228  shape[3] = tau3_num/tau_den;
1229  shape[4] = dR2hardest;
1230  shape[5] = dPhi2hardest;
1231  }
1232 
1233  shape[6] = 10.;
1234  return shape;
1235 
1236  } while(0);
1237 
1238  shape[6] = -10.;
1239  return shape;
1240 }
1241 
AliAnalysisTaskEmcalMissingEnergy::N_subjettiness
Float_t * N_subjettiness(AliEmcalJet *mainJet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1159
AliAnalysisTaskEmcalMissingEnergy::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:950
AliEmcalJetFinder::SetRadius
void SetRadius(Double_t val)
Definition: AliEmcalJetFinder.h:40
AliAnalysisTaskEmcalMissingEnergy::Circularity
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:642
AliAnalysisTaskEmcalMissingEnergy::GetJetMass
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:563
AliPicoTrack.h
Double_t
double Double_t
Definition: External.C:58
TH3F
Definition: External.C:260
AliEmcalJetFinder::SetJetMinPt
void SetJetMinPt(Double_t val)
Definition: AliEmcalJetFinder.h:45
AliEmcalJetFinder.h
AliEmcalJetShapeProperties::GetSecondOrderSubtractedSigma2
Double_t GetSecondOrderSubtractedSigma2() const
Definition: AliEmcalJetShapeProperties.h:93
TH2F
Definition: External.C:236
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:377
AliEmcalJetShapeProperties::GetSecondOrderSubtractedConstituent
Double_t GetSecondOrderSubtractedConstituent() const
Definition: AliEmcalJetShapeProperties.h:104
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:121
AliEmcalJet::Py
Double_t Py() const
Definition: AliEmcalJet.h:107
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:117
AliEmcalPythiaInfo::GetPythiaEventWeight
Float_t GetPythiaEventWeight() const
Definition: AliEmcalPythiaInfo.h:56
AliAnalysisTaskEmcalMissingEnergy::RelativePhiFancy
Double_t RelativePhiFancy(Double_t mphi, Double_t vphi)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1017
AliEmcalJet::GetLabel
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
AliAnalysisTaskEmcalMissingEnergy::Tau2Num
Float_t Tau2Num(AliEmcalJet *jet, AliEmcalJet *subJet1hardest, AliEmcalJet *subJet2hardest, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1121
AliEmcalParticle.h
AliAnalysisTaskEmcalMissingEnergy::Run
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.
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:329
AliAnalysisTaskEmcalMissingEnergy::GetJetAngularity
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:598
AliAnalysisTaskEmcalMissingEnergy::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:337
AliAnalysisTaskEmcalMissingEnergy::GetJetpTD
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:632
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
AliEmcalJet::TrackAt
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
AliAnalysisTaskEmcalMissingEnergy::Tau1Num
Float_t Tau1Num(AliEmcalJet *jet, AliEmcalJet *subJet1hardest, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1081
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:22
AliEmcalJet::Px
Double_t Px() const
Definition: AliEmcalJet.h:106
AliAnalysisTaskEmcalMissingEnergy::LeSub
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:730
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1623
AliJetContainer::GetParticleContainer
AliParticleContainer * GetParticleContainer() const
Definition: AliJetContainer.h:180
AliPicoTrack::Eta
Double_t Eta() const
Definition: AliPicoTrack.h:37
AliAnalysisTaskEmcalMissingEnergy::GetJetLeSub
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:767
AliEmcalJetFinder::SetJetAlgorithm
void SetJetAlgorithm(Int_t val)
Definition: AliEmcalJetFinder.h:41
AliAnalysisTaskEmcalMissingEnergy::GetJetCircularity
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:720
Int_t
int Int_t
Definition: External.C:63
AliEmcalParticle::Eta
Double_t Eta() const
Definition: AliEmcalParticle.h:36
AliEmcalParticle::Pt
Double_t Pt() const
Definition: AliEmcalParticle.h:24
UInt_t
unsigned int UInt_t
Definition: External.C:33
Float_t
float Float_t
Definition: External.C:68
AliParticleContainer.h
AliJetContainer::GetNJets
Int_t GetNJets() const
Definition: AliJetContainer.h:153
AliEmcalJetShapeProperties::GetSecondOrderSubtractedAngularity
Double_t GetSecondOrderSubtractedAngularity() const
Definition: AliEmcalJetShapeProperties.h:63
AliJetContainer.h
AliPicoTrack::Phi
Double_t Phi() const
Definition: AliPicoTrack.h:33
AliAnalysisTaskEmcal::fCent
Double_t fCent
!event centrality
Definition: AliAnalysisTaskEmcal.h:886
AliRhoParameter.h
AliPicoTrack::GetTrackType
Byte_t GetTrackType() const
Definition: AliPicoTrack.h:43
AliPicoTrack::Pt
Double_t Pt() const
Definition: AliPicoTrack.h:23
AliAnalysisTaskEmcalMissingEnergy::Tau3Num
Float_t Tau3Num(AliEmcalJet *jet, AliEmcalJet *subJet1hardest, AliEmcalJet *subJet2hardest, AliEmcalJet *subJet3hardest, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1139
AliEmcalJet::GetZ
Double_t GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz) const
Definition: AliEmcalJet.cxx:458
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:299
Data
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)
Definition: charmCutsOptimization.C:371
AliEmcalJetFinder::GetJet
AliEmcalJet * GetJet(Int_t index)
Definition: AliEmcalJetFinder.h:37
AliAnalysisTaskEmcalMissingEnergy::R_distance
Double_t R_distance(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:999
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliAnalysisTaskEmcalMissingEnergy::TauDen
Float_t TauDen(AliEmcalJet *mainJet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1067
AliAnalysisTaskEmcalMissingEnergy::GetSigma2
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:851
AliAnalysisTaskEmcalMissingEnergy::RelativePhi
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:936
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:909
AliAnalysisTaskEmcalMissingEnergy::GetJetNumberOfConstituents
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:778
AliEmcalJetShapeProperties::GetSecondOrderSubtractedCircularity
Double_t GetSecondOrderSubtractedCircularity() const
Definition: AliEmcalJetShapeProperties.h:83
AliAnalysisTaskEmcalMissingEnergy::SelectTrigger
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:862
AliEmcalPythiaInfo
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
Definition: AliEmcalPythiaInfo.h:31
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:256
AliAnalysisTaskEmcalMissingEnergy::Tau1Num_full
Float_t Tau1Num_full(AliEmcalJet *mainJet, AliEmcalJetFinder *finder, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:1097
AliAnalysisTaskEmcalJet
Base task in the EMCAL jet framework.
Definition: AliAnalysisTaskEmcalJet.h:30
AliEmcalJet
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
AliAnalysisTaskEmcal::GetPythiaInfo
const AliEmcalPythiaInfo * GetPythiaInfo() const
Definition: AliAnalysisTaskEmcal.h:333
AliEmcalJet::Pz
Double_t Pz() const
Definition: AliEmcalJet.h:108
AliJetContainer::GetFractionSharedPt
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Definition: AliJetContainer.cxx:806
AliEmcalPythiaInfo.h
Declaration of class AliEmcalPythiaInfo.
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:336
Bool_t
bool Bool_t
Definition: External.C:53
AliEmcalJet::GetShapeProperties
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
AliAnalysisTaskEmcalMissingEnergy::Angularity
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:572
AliAnalysisTaskEmcalMissingEnergy::Sigma2
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:790
AliEmcalJet::M
Double_t M() const
Definition: AliEmcalJet.h:120
AliAnalysisTaskEmcalMissingEnergy::PTD
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
Definition: AliAnalysisTaskEmcalMissingEnergy.cxx:609
AliEmcalJet.h
AliJetContainer
Container for jet within the EMCAL jet framework.
Definition: AliJetContainer.h:37
AliEmcalJet::GetPtSortedTrackConstituentIndexes
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Definition: AliEmcalJet.cxx:424
AliJetContainer::GetJet
AliEmcalJet * GetJet(Int_t i) const
Definition: AliJetContainer.cxx:270
AliAnalysisTaskEmcalMissingEnergy.h