AliPhysics  5364b50 (5364b50)
AliAnalysisTaskRecursiveSoftDrop.cxx
Go to the documentation of this file.
1 //
2 // Basic analysis task.
3 //
4 // Basic analysis task template for analysis jets storing information in both tree
5 // branches and histograms
6 
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <TCanvas.h>
12 #include <THnSparse.h>
13 #include <TTree.h>
14 #include <TList.h>
15 #include <TLorentzVector.h>
16 #include <TProfile.h>
17 #include <TChain.h>
18 #include <TSystem.h>
19 #include <TFile.h>
20 #include <TKey.h>
21 #include <AliAnalysisDataSlot.h>
22 #include <AliAnalysisDataContainer.h>
23 #include <vector>
24 #include "TMatrixD.h"
25 #include "TMatrixDSym.h"
26 #include "TMatrixDSymEigen.h"
27 #include "TVector3.h"
28 #include "TVector2.h"
29 #include "AliVCluster.h"
30 #include "AliVTrack.h"
31 #include "AliEmcalJet.h"
32 #include "AliRhoParameter.h"
33 #include "AliLog.h"
34 #include "AliEmcalParticle.h"
35 #include "AliMCEvent.h"
36 #include "AliGenPythiaEventHeader.h"
37 #include "AliAODMCHeader.h"
38 #include "AliMCEvent.h"
39 #include "AliAnalysisManager.h"
40 #include "AliJetContainer.h"
41 #include "AliParticleContainer.h"
42 //#include "AliPythiaInfo.h"
43 #include "TRandom3.h"
44 #include "AliPicoTrack.h"
45 #include "AliEmcalJetFinder.h"
46 #include "AliAODEvent.h"
48 
49 #include "FJ_includes.h"
50 
51 //Globals
52 using std::cout;
53 using std::endl;
54 
56 
57 //________________________________________________________________________
60  fContainer(0),
61  fJetShapeSub(kNoSub),
62  fPtThreshold(-9999.),
63  fCentSelectOn(kTRUE),
64  fCentMin(0),
65  fCentMax(10),
66  fJetRadius(0.4),
67  fSharedFractionPtMin(0.5),
68  fhJetPt(0x0),
69  fhJetPhi(0x0),
70  fhJetEta(0x0),
71  fhDetJetPt_Matched(0x0),
72  fTreeRecursive_Det(0),
73  fTreeRecursive_True(0)
74 
75 {
76  for(Int_t i=0;i<4;i++){
77  fShapesVar_Det[i]=0;
78  fShapesVar_True[i]=0;
79  }
80  SetMakeGeneralHistograms(kTRUE);
81  DefineOutput(1, TList::Class());
82  DefineOutput(2, TTree::Class());
83  DefineOutput(3, TTree::Class());
84 }
85 
86 //________________________________________________________________________
88  AliAnalysisTaskEmcalJet(name, kTRUE),
89  fContainer(0),
90  fJetShapeSub(kNoSub),
91  fPtThreshold(-9999.),
92  fCentSelectOn(kTRUE),
93  fCentMin(0),
94  fCentMax(10),
95  fJetRadius(0.4),
96  fSharedFractionPtMin(0.5),
97  fhJetPt(0x0),
98  fhJetPhi(0x0),
99  fhJetEta(0x0),
100  fhDetJetPt_Matched(0x0),
101  fTreeRecursive_Det(0),
102  fTreeRecursive_True(0)
103 {
104  // Standard constructor.
105  for(Int_t i=0;i<4;i++){
106  fShapesVar_Det[i]=0;
107  fShapesVar_True[i]=0;
108  }
110  DefineOutput(1, TList::Class());
111  DefineOutput(2, TTree::Class());
112  DefineOutput(3, TTree::Class());
113 }
114 
115 //________________________________________________________________________
117 {
118  // Destructor.
119 }
120 
121 //________________________________________________________________________
123 {
124  // Create user output.
125 
127 
128  Bool_t oldStatus = TH1::AddDirectoryStatus();
129  TH1::AddDirectory(kFALSE);
130  TH1::AddDirectory(oldStatus);
131  //create a tree used for the MC data and making a 4D response matrix
132  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
133  fTreeRecursive_Det = new TTree(nameoutput, nameoutput);
134  const char* nameoutput2 = GetOutputSlot(3)->GetContainer()->GetName();
135  fTreeRecursive_True = new TTree(nameoutput2, nameoutput2);
136 
137  const Int_t intBranches = 4;
138 
139  std::vector<TString> fShapesVarNames_Det(intBranches), fShapesVarNames_True(intBranches);
140 
141 
142  fShapesVarNames_Det[0] = "Pt";
143  fShapesVarNames_Det[1] = "Z";
144  fShapesVarNames_Det[2] = "Theta";
145  fShapesVarNames_Det[3] = "N";
146  fShapesVarNames_True[0] = "Pt_Truth";
147  fShapesVarNames_True[1] = "Z_Truth";
148  fShapesVarNames_True[2] = "Theta_Truth";
149  fShapesVarNames_True[3] = "N_Truth";
150 
151  for(Int_t ivar=0; ivar < intBranches; ivar++){
152  cout<<"looping over variables"<<endl;
153  fTreeRecursive_Det->Branch(fShapesVarNames_Det[ivar].Data(), &fShapesVar_Det[ivar], Form("%s/D", fShapesVarNames_Det[ivar].Data()));
154  fTreeRecursive_True->Branch(fShapesVarNames_True[ivar].Data(), &fShapesVar_True[ivar], Form("%s/D", fShapesVarNames_True[ivar].Data()));
155 
156  }
157 
158  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
159  fOutput->Add(fhJetPt);
160  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
161  fOutput->Add(fhJetPhi);
162  fhJetEta= new TH1F("fhJetEta", "Jet Eta",100,-2,2);
163  fOutput->Add(fhJetEta);
164  fhDetJetPt_Matched= new TH1F("fhDetJetPt_Matched", "Jet Pt",200,-0.5,199.5 );
166 
167  PostData(1,fOutput);
168  PostData(2,fTreeRecursive_Det);
169  PostData(3,fTreeRecursive_True);
170  // delete [] fShapesVarNames;
171 }
172 
173 //________________________________________________________________________
175 {
176  // Run analysis code here, if needed. It will be executed before FillHistograms().
177 
178 
179  return kTRUE;
180 }
181 
182 //________________________________________________________________________
184 {
185  if (fCentSelectOn){
186  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
187  }
188 
189  if(fJetType == kData){
190  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
191  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
192  Double_t JetPhi=0;
193  Double_t JetPt_ForThreshold=0;
194  if(JetCont) {
195  JetCont->ResetCurrentID();
196  while((Jet1=JetCont->GetNextAcceptJet())) {
197  if(!Jet1) continue;
198  if (fJetShapeSub==kNoSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
199  else JetPt_ForThreshold = Jet1->Pt();
200  if(JetPt_ForThreshold<fPtThreshold) {
201  continue;
202  }
203  else {
204  fhJetPt->Fill(Jet1->Pt());
205  JetPhi=Jet1->Phi();
206  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
207  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
208  fhJetPhi->Fill(JetPhi);
209  fhJetEta->Fill(Jet1->Eta());
210  RecursiveParents(Jet1,JetCont,1,kFALSE); //Third argument = reclustering algorithm (0=Antikt,1=CA,2=kt)
211  }
212  }
213  }
214  }
215 
216  if(fJetType == kEmb){
217  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
218  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet //For matching SubtractedHybrid->DetPythia this jet container is also Subtracted Hybrid
219  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
220  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
221  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
222  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
223  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
224  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
225 
226 
227 
228  Bool_t JetsMatched = kFALSE;
229  Double_t JetPtThreshold;
230  JetContHybridS->ResetCurrentID();
231  JetContHybridUS->ResetCurrentID();
232  JetContPythDet->ResetCurrentID();
233  JetContPythTrue->ResetCurrentID();
234 
235  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){ //Get next Subtracted hybrid jet
236  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
237  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
238  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue; //check pT is above threshold
239  Int_t JetNumber=-1;
240  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
241  JetHybridUS = JetContHybridUS->GetJet(i); //Get unsubtracted jets in order
242  if (!JetHybridUS) continue;
243 
244  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) { //check if it mataches with current subtracted hybrid jet
245  JetNumber=i;
246  }
247  }
248  if(JetNumber==-1) continue;
249  JetHybridUS=JetContHybridUS->GetJet(JetNumber); //Get the matched Unsubtracted jet
250  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin) { //Check that the US closest jet shares a minimum
251  continue; // pT with Detector level pythia jet
252  }
253  JetPythDet=JetHybridUS->ClosestJet(); //get the closest jet that has passed the shared pT cut
254  if (!JetHybridUS) {
255  Printf("Unsubtracted embedded jet does not exist, returning");
256  continue;
257  }
258  if (!JetPythDet) continue;
259  UInt_t rejectionReason = 0;
260  if (!(JetContPythDet->AcceptJet(JetPythDet,rejectionReason))) continue; //Check the detector level just is accepted
261  fhDetJetPt_Matched->Fill(JetPythDet->Pt()); //Fill only matched detector level jets for tagging efficiency comparison
262  JetPythTrue=JetPythDet->ClosestJet(); //Get the corresponding pythia true jet
263  if(!JetPythTrue) continue;
264  JetsMatched=kTRUE; //jets have been matched
265 
266 
267  fhJetPt->Fill(JetHybridS->Pt());
268  Double_t JetPhi=JetHybridS->Phi();
269  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
270  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
271  fhJetPhi->Fill(JetPhi);
272  fhJetEta->Fill(JetHybridS->Eta());
273  RecursiveParents(JetHybridS,JetContHybridS,1,kFALSE); //Third argument = reclustering algorithm (0=Antikt,1=CA,2=kt)
274  RecursiveParents(JetPythTrue,JetContPythTrue,1,kTRUE); //Third argument = reclustering algorithm (0=Antikt,1=CA,2=kt)
275 
276 
277 
278  }
279  }
280 
281 
282  return kTRUE;
283 }
284 
285 //_________________________________________________________________________
287  std::vector<fastjet::PseudoJet> fInputVectors;
288  fInputVectors.clear();
289  fastjet::PseudoJet PseudoTracks;
290  double xflagalgo=0;
291  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
292 
293  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
294  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
295  if (!fTrk) continue;
296  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
297  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
298  fInputVectors.push_back(PseudoTracks);
299 
300  }
301 
302 
303 
304  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
305  if(ReclusterAlgo==0){ xflagalgo=0.5;
306  jetalgo=fastjet::kt_algorithm ;}
307 
308  if(ReclusterAlgo==1){ xflagalgo=1.5;
309  jetalgo=fastjet::cambridge_algorithm;
310  }
311  if(ReclusterAlgo==2){ xflagalgo=2.5;
312  jetalgo=fastjet::antikt_algorithm;
313  }
314 
315  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
316 
317  try {
318  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
319  std::vector<fastjet::PseudoJet> fOutputJets;
320  fOutputJets.clear();
321  fOutputJets=fClustSeqSA.inclusive_jets(0);
322 
323  fastjet::PseudoJet jj;
324  fastjet::PseudoJet j1;
325  fastjet::PseudoJet j2;
326  jj=fOutputJets[0];
327  Int_t n = 0;
328  Double_t jet_pT=jj.perp();
329  while(jj.has_parents(j1,j2)){
330  n++;
331  if(j1.perp() < j2.perp()) swap(j1,j2);
332  double delta_R=j1.delta_R(j2);
333  double z=j2.perp()/(j1.perp()+j2.perp());
334  if(bTruth) {
335  fShapesVar_True[0]=jet_pT;
336  fShapesVar_True[1]=z;
337  fShapesVar_True[2]=delta_R;
338  fShapesVar_True[3]=n;
339  fTreeRecursive_True->Fill();
340  }
341  else {
342  fShapesVar_Det[0]=jet_pT;
343  fShapesVar_Det[1]=z;
344  fShapesVar_Det[2]=delta_R;
345  fShapesVar_Det[3]=n;
346  fTreeRecursive_Det->Fill();
347  }
348  jj=j1;
349  }
350 
351  } catch (fastjet::Error) {
352  AliError(" [w] FJ Exception caught.");
353  //return -1;
354  }
355  return;
356 }
357 
358 
359 //________________________________________________________________________
361  //
362  // retrieve event objects
363  //
365  return kFALSE;
366 
367  return kTRUE;
368 }
369 
370 
371 //_______________________________________________________________________
373 {
374  // Called once at the end of the analysis.
375 
376 }
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Double_t Area() const
Definition: AliEmcalJet.h:130
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo, Bool_t bTruth)
double Double_t
Definition: External.C:58
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Bool_t FillHistograms()
Function filling histograms.
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
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.
Int_t GetNJets() const
Double_t fCent
!event centrality
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)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Double_t GetRhoVal(Int_t i=0) const
AliEmcalList * fOutput
!output list
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void swap(AliEmcalContainerIndexMap< X, Y > &first, AliEmcalContainerIndexMap< X, Y > &second)
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
Container for jet within the EMCAL jet framework.
AliEmcalJet * GetJet(Int_t i) const