AliPhysics  58f3d52 (58f3d52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  fhJetPt(0x0),
68  fhJetPhi(0x0),
69  fhJetEta(0x0),
70  fTreeRecursive_Det(0),
71  fTreeRecursive_True(0)
72 
73 {
74  for(Int_t i=0;i<4;i++){
75  fShapesVar_Det[i]=0;
76  fShapesVar_True[i]=0;
77  }
78  SetMakeGeneralHistograms(kTRUE);
79  DefineOutput(1, TList::Class());
80  DefineOutput(2, TTree::Class());
81  DefineOutput(3, TTree::Class());
82 }
83 
84 //________________________________________________________________________
86  AliAnalysisTaskEmcalJet(name, kTRUE),
87  fContainer(0),
88  fJetShapeSub(kNoSub),
89  fPtThreshold(-9999.),
90  fCentSelectOn(kTRUE),
91  fCentMin(0),
92  fCentMax(10),
93  fJetRadius(0.4),
94  fhJetPt(0x0),
95  fhJetPhi(0x0),
96  fhJetEta(0x0),
97  fTreeRecursive_Det(0),
98  fTreeRecursive_True(0)
99 {
100  // Standard constructor.
101  for(Int_t i=0;i<4;i++){
102  fShapesVar_Det[i]=0;
103  fShapesVar_True[i]=0;
104  }
106  DefineOutput(1, TList::Class());
107  DefineOutput(2, TTree::Class());
108  DefineOutput(3, TTree::Class());
109 }
110 
111 //________________________________________________________________________
113 {
114  // Destructor.
115 }
116 
117 //________________________________________________________________________
119 {
120  // Create user output.
121 
123 
124  Bool_t oldStatus = TH1::AddDirectoryStatus();
125  TH1::AddDirectory(kFALSE);
126  TH1::AddDirectory(oldStatus);
127  //create a tree used for the MC data and making a 4D response matrix
128  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
129  fTreeRecursive_Det = new TTree(nameoutput, nameoutput);
130  const char* nameoutput2 = GetOutputSlot(3)->GetContainer()->GetName();
131  fTreeRecursive_True = new TTree(nameoutput2, nameoutput2);
132 
133  const Int_t intBranches = 4;
134 
135  std::vector<TString> fShapesVarNames_Det(intBranches), fShapesVarNames_True(intBranches);
136 
137 
138  fShapesVarNames_Det[0] = "Pt";
139  fShapesVarNames_Det[1] = "Z";
140  fShapesVarNames_Det[2] = "Theta";
141  fShapesVarNames_Det[3] = "N";
142  fShapesVarNames_True[0] = "Pt_Truth";
143  fShapesVarNames_True[1] = "Z_Truth";
144  fShapesVarNames_True[2] = "Theta_Truth";
145  fShapesVarNames_True[3] = "N_Truth";
146 
147  for(Int_t ivar=0; ivar < intBranches; ivar++){
148  cout<<"looping over variables"<<endl;
149  fTreeRecursive_Det->Branch(fShapesVarNames_Det[ivar].Data(), &fShapesVar_Det[ivar], Form("%s/D", fShapesVarNames_Det[ivar].Data()));
150  fTreeRecursive_True->Branch(fShapesVarNames_True[ivar].Data(), &fShapesVar_True[ivar], Form("%s/D", fShapesVarNames_True[ivar].Data()));
151 
152  }
153 
154  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
155  fOutput->Add(fhJetPt);
156  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
157  fOutput->Add(fhJetPhi);
158  fhJetEta= new TH1F("fhJetEta", "Jet Eta",100,-2,2);
159  fOutput->Add(fhJetEta);
160 
161  PostData(1,fOutput);
162  PostData(2,fTreeRecursive_Det);
163  PostData(3,fTreeRecursive_True);
164  // delete [] fShapesVarNames;
165 }
166 
167 //________________________________________________________________________
169 {
170  // Run analysis code here, if needed. It will be executed before FillHistograms().
171 
172 
173  return kTRUE;
174 }
175 
176 //________________________________________________________________________
178 {
179  if (fCentSelectOn){
180  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
181  }
182  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
183  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
184  Double_t JetPhi=0;
185  Double_t JetPt_ForThreshold=0;
186  if(JetCont) {
187  JetCont->ResetCurrentID();
188  while((Jet1=JetCont->GetNextAcceptJet())) {
189  if(!Jet1) continue;
190  if (fJetShapeSub==kNoSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
191  else JetPt_ForThreshold = Jet1->Pt();
192  if(JetPt_ForThreshold<fPtThreshold) {
193  continue;
194  }
195  else {
196  fhJetPt->Fill(Jet1->Pt());
197  JetPhi=Jet1->Phi();
198  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
199  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
200  fhJetPhi->Fill(JetPhi);
201  fhJetEta->Fill(Jet1->Eta());
202  RecursiveParents(Jet1,JetCont,1,kFALSE); //Third argument = reclustering algorithm (0=Antikt,1=CA,2=kt)
203  }
204  }
205  }
206 
207  return kTRUE;
208 }
209 
210 //_________________________________________________________________________
212  std::vector<fastjet::PseudoJet> fInputVectors;
213  fInputVectors.clear();
214  fastjet::PseudoJet PseudoTracks;
215  double xflagalgo=0;
216  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
217 
218  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
219  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
220  if (!fTrk) continue;
221  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
222  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
223  fInputVectors.push_back(PseudoTracks);
224 
225  }
226 
227 
228 
229  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
230  if(ReclusterAlgo==0){ xflagalgo=0.5;
231  jetalgo=fastjet::kt_algorithm ;}
232 
233  if(ReclusterAlgo==1){ xflagalgo=1.5;
234  jetalgo=fastjet::cambridge_algorithm;
235  }
236  if(ReclusterAlgo==2){ xflagalgo=2.5;
237  jetalgo=fastjet::antikt_algorithm;
238  }
239 
240  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
241 
242  try {
243  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
244  std::vector<fastjet::PseudoJet> fOutputJets;
245  fOutputJets.clear();
246  fOutputJets=fClustSeqSA.inclusive_jets(0);
247 
248  fastjet::PseudoJet jj;
249  fastjet::PseudoJet j1;
250  fastjet::PseudoJet j2;
251  jj=fOutputJets[0];
252  Int_t n = 0;
253  Double_t jet_pT=jj.perp();
254  while(jj.has_parents(j1,j2)){
255  n++;
256  if(j1.perp() < j2.perp()) swap(j1,j2);
257  double delta_R=j1.delta_R(j2);
258  double z=j2.perp()/(j1.perp()+j2.perp());
259  if(bTruth) {
260  fShapesVar_True[0]=jet_pT;
261  fShapesVar_True[1]=z;
262  fShapesVar_True[2]=delta_R;
263  fShapesVar_True[3]=n;
264  fTreeRecursive_True->Fill();
265  }
266  else {
267  fShapesVar_Det[0]=jet_pT;
268  fShapesVar_Det[1]=z;
269  fShapesVar_Det[2]=delta_R;
270  fShapesVar_Det[3]=n;
271  fTreeRecursive_Det->Fill();
272  }
273  jj=j1;
274  }
275 
276  } catch (fastjet::Error) {
277  AliError(" [w] FJ Exception caught.");
278  //return -1;
279  }
280  return;
281 }
282 
283 
284 //________________________________________________________________________
286  //
287  // retrieve event objects
288  //
290  return kFALSE;
291 
292  return kTRUE;
293 }
294 
295 
296 //_______________________________________________________________________
298 {
299  // Called once at the end of the analysis.
300 
301 }
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
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
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
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 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.