AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetHadCorQA.cxx
Go to the documentation of this file.
2 
3 #include <TCanvas.h>
4 #include <TChain.h>
5 #include <TClonesArray.h>
6 #include <TH1F.h>
7 #include <TH2F.h>
8 #include <TH3F.h>
9 #include <THnSparse.h>
10 #include <TList.h>
11 #include <TLorentzVector.h>
12 #include <TParameter.h>
13 #include <TParticle.h>
14 #include <TTree.h>
15 #include <TVector3.h>
16 
17 #include "AliAODEvent.h"
18 #include "AliAnalysisManager.h"
19 #include "AliAnalysisTask.h"
20 #include "AliCentrality.h"
21 #include "AliESDEvent.h"
22 #include "AliESDInputHandler.h"
23 #include "AliEmcalJet.h"
24 #include "AliVCluster.h"
25 #include "AliVTrack.h"
26 #include "AliRhoParameter.h"
27 #include "AliEmcalParticle.h"
28 #include "AliPicoTrack.h"
29 #include "AliEMCALGeometry.h"
30 using std::vector;
31 
33 
34 //________________________________________________________________________
36  AliAnalysisTaskEmcalJet("jethadcor",kFALSE),
37  fCalo2Name(),
38  fCaloClusters2(),
39  fMCParticlesName(),
40  fMCParticles(),
41  fHistRhovsCent(0),
42  fHistNjetvsCent(0),
43  fHistNTMatchvsPtvsNtack0(0)
44 {
45  for (int i = 0;i<3;i++){
46  fHistNEFvsPt[i] = 0;
47  fHistNTMatchvsPt[i] = 0;
48  fHistNCMatchvsPt[i] = 0;
49  fHistHadCorvsPt[i] = 0;
50  fHistNEFvsPtBias[i] = 0;
51  fHistNconvsPt[i] = 0;
52  fHistNtvsPt[i] = 0;
53  fHistNcvsPt[i] = 0;
54  fHistNTMatchvsPtBias[i] = 0;
55  fHistNCMatchvsPtBias[i] = 0;
56  fHistHadCorvsPtBias[i] = 0;
57  fHistNconvsPtBias[i] = 0;
58  fHistNtvsPtBias[i] = 0;
59  fHistNcvsPtBias[i] = 0;
60  }
61  // Default constructor.
62 
63  SetMakeGeneralHistograms(kTRUE);
64 }
65 
66 //________________________________________________________________________
68  AliAnalysisTaskEmcalJet(name,kTRUE),
69  fCalo2Name(),
70  fCaloClusters2(),
71  fMCParticlesName(),
72  fMCParticles(),
73  fHistRhovsCent(0),
74  fHistNjetvsCent(0),
75  fHistNTMatchvsPtvsNtack0(0)
76  {
77  for (int i = 0;i<3;i++){
78  fHistNEFvsPt[i] = 0;
79  fHistNTMatchvsPt[i] = 0;
80  fHistNCMatchvsPt[i] = 0;
81  fHistHadCorvsPt[i] = 0;
82  fHistNEFvsPtBias[i] = 0;
83  fHistNconvsPt[i] = 0;
84  fHistNtvsPt[i] = 0;
85  fHistNcvsPt[i] = 0;
86  fHistNTMatchvsPtBias[i] = 0;
87  fHistNCMatchvsPtBias[i] = 0;
88  fHistHadCorvsPtBias[i] = 0;
89  fHistNconvsPtBias[i] = 0;
90  fHistNtvsPtBias[i] = 0;
91  fHistNcvsPtBias[i] = 0;
92  }
94  }
95 
96 //________________________________________________________________________
98 {
99  if (! fCreateHisto)
100  return;
102  fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 500, 0, 500);
103  fHistNjetvsCent = new TH2F("NjetvsCent", "NjetvsCent", 100, 0.0, 100.0, 100, 0, 100);
104 
105  for (int i = 0;i<3;i++){
106  char name[200];
107  TString title;
108  sprintf(name,"NEFvsPt%i",i);
109  fHistNEFvsPt[i] = new TH2F(name, name, 100,0,1,500,-250,250);
110  fOutput->Add(fHistNEFvsPt[i]);
111  sprintf(name,"NTMatchvsPt%i",i);
112  fHistNTMatchvsPt[i] = new TH2F(name, name, 100,0,100,500,-250,250);
113  fOutput->Add(fHistNTMatchvsPt[i]);
114  sprintf(name,"NCMatchvsPt%i",i);
115  fHistNCMatchvsPt[i] = new TH2F(name, name, 100,0,100,500,-250,250);
116  fOutput->Add(fHistNCMatchvsPt[i]);
117  sprintf(name,"HadCorvsPt%i",i);
118  fHistHadCorvsPt[i] = new TH2F(name, name, 1000,0,500,500,-250,250);
119  fOutput->Add(fHistHadCorvsPt[i]);
120  sprintf(name,"NconvsPt%i",i);
121  fHistNconvsPt[i] = new TH2F(name, name, 200,0,200,500,-250,250);
122  fOutput->Add(fHistNconvsPt[i]);
123  sprintf(name,"NtvsPt%i",i);
124  fHistNtvsPt[i] = new TH2F(name, name, 200,0,200,500,-250,250);
125  fOutput->Add(fHistNtvsPt[i]);
126  sprintf(name,"NcvsPt%i",i);
127  fHistNcvsPt[i] = new TH2F(name, name, 200,0,200,500,-250,250);
128  fOutput->Add(fHistNcvsPt[i]);
129  sprintf(name,"NEFvsPtBias%i",i);
130  fHistNEFvsPtBias[i] = new TH2F(name, name, 100,0,1,500,-250,250);
131  fOutput->Add(fHistNEFvsPtBias[i]);
132  sprintf(name,"NTMatchvsPtBias%i",i);
133  fHistNTMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
134  fOutput->Add(fHistNTMatchvsPtBias[i]);
135  sprintf(name,"NCMatchvsPtBias%i",i);
136  fHistNCMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
137  fOutput->Add(fHistNCMatchvsPtBias[i]);
138  sprintf(name,"HadCorvsPtBias%i",i);
139  fHistHadCorvsPtBias[i] = new TH2F(name, name, 1000,0,500,500,-250,250);
140  fOutput->Add(fHistHadCorvsPtBias[i]);
141  sprintf(name,"NconvsPtBias%i",i);
142  fHistNconvsPtBias[i] = new TH2F(name, name, 200,0,200,500,-250,250);
143  fOutput->Add(fHistNconvsPtBias[i]);
144  sprintf(name,"NtvsPtBias%i",i);
145  fHistNtvsPtBias[i] = new TH2F(name, name, 200,0,200,500,-250,250);
146  fOutput->Add(fHistNtvsPtBias[i]);
147  sprintf(name,"NcvsPtBias%i",i);
148  fHistNcvsPtBias[i] = new TH2F(name, name, 200,0,200,500,-250,250);
149  fOutput->Add(fHistNcvsPtBias[i]);
150  }
151  fHistNTMatchvsPtvsNtack0 = new TH3F("NTMmatchvsPtvsNtrack0", "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);
152 
153  // for (Int_t i = 0;i<6;++i){
154  // name = TString(Form("JetPtvsTrackPt_%i",i));
155  // title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
156  // fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
157  // fOutput->Add(fHistJetPtvsTrackPt[i]);
158 
159  // }
160  fOutput->Add(fHistRhovsCent);
161  fOutput->Add(fHistNjetvsCent);
163  PostData(1, fOutput);
164 }
165 
166 //________________________________________________________________________
167 
169 {
170  // Get centrality bin.
171 
172  Int_t centbin = -1;
173  if (cent>=0 && cent<10)
174  centbin = 0;
175  else if (cent>=10 && cent<30)
176  centbin = 1;
177  else if (cent>=30 && cent<50)
178  centbin = 2;
179  else if (cent>50)
180  centbin =3;
181  return centbin;
182 }
183 
184 //________________________________________________________________________
185 
187 {
188  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
189  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
190  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
191  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
192  double dphi = mphi-vphi;
193  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
194  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
195 
196  return dphi;//dphi in [-Pi, Pi]
197 }
198 
199 //________________________________________________________________________
201 {
203 
204 // AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
205 // if (fCaloName == "CaloClusters")
206 // am->LoadBranch("CaloClusters");
207 
208  if (!fCalo2Name.IsNull() && !fCaloClusters2){
209  fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
210  if (!fCaloClusters2){
211  AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
212  fLocalInitialized = kFALSE;
213  return;
214  }
215  }
216 
217  if (!fMCParticlesName.IsNull() && !fMCParticles){
218  fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
219  if (!fMCParticles){
220  AliError(Form("%s: Could not retrieve MC Particles %s!",GetName(),fMCParticlesName.Data()));
221  fLocalInitialized = kFALSE;
222  return;
223  }
224  }
225  //TString fMCParticlesName;
226  // TClonesArray *fMCParticles;
227  // fCaloClusters2(),
228  // else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
229 // AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
230 // fCaloClusters2 = 0;
231 // fLocalInitialized = kFALSE;
232 // return;
233  // }
234 }
235 
236 
237 //________________________________________________________________________
239 {
240  Int_t centbin = GetCentBin(fCent);
241  //for pp analyses we will just use the first centrality bin
242  if (GetBeamType()==0)
243  centbin = 0;
244  if (centbin>2)
245  return kTRUE;
246  if (!fTracks)
247  return kTRUE;
248  if (!fCaloClusters)
249  return kTRUE;
250  if (!fCaloClusters2)
251  return kTRUE;
252  const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
253  const Int_t nTrack = fTracks->GetEntriesFast();
254  // const Int_t nMC = fMCParticles->GetEntriesFast();
255 
256  TString fRhoScaledName = fRhoName;
257  fRho = GetRhoFromEvent(fRhoScaledName);
258  fRhoVal = fRho->GetVal();
259  if (GetBeamType()==0)
260  fRhoVal = 0;
262  const Int_t Njets = fJets->GetEntriesFast();
263 
264  Int_t NjetAcc = 0;
265 
266 
267  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
268  Int_t TrackMatch = 0;
269  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
270  if (!jet)
271  continue;
272  if (jet->Area()==0)
273  continue;
274  if (jet->Pt()<0.1)
275  continue;
276  if (jet->MaxTrackPt()>100)
277  continue;
278  if (! AcceptJet(jet))
279  continue;
280  NjetAcc++;
281  vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
282  vector<Int_t> cluster_id2;
283  Double_t Esub = 0; //total E subtracted from the jet
284  Double_t jetPt = -500;
285  jetPt = jet->Pt()-jet->Area()*fRhoVal;
286  for (int i = 0;i<nTrack;i++){
287  AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(i));
288  if (!emctrack)
289  continue;
290  if (!emctrack->IsEMCAL())
291  continue;
292  AliVTrack *track = (AliVTrack*)emctrack->GetTrack();
293  if (! track)
294  continue;
295  if (! AcceptTrack(track))
296  continue;
297  if (! IsJetTrack(jet,i,false))
298  continue;
299  Int_t iClus = track->GetEMCALcluster();
300  if (iClus<0)
301  continue;
302  //we have the id of the matched cluster of a track constituent
303  // cout<<"track label for matched,accepted jet track is "<<track->GetLabel()<<endl;
304  bool ischecked = false;
305  for (UInt_t icid = 0;icid<cluster_id.size();icid++)
306  if (cluster_id[icid] == iClus)
307  ischecked = true; // we've already looked at this uncorrected cluster
308  if (ischecked)
309  continue; //no need to go further
310  AliEmcalParticle *emcluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
311  AliVCluster *cluster = emcluster->GetCluster();
312  if (! cluster)
313  continue;
314  if (! AcceptCluster(cluster))
315  continue;
316  Double_t etadiff = 999;
317  Double_t phidiff = 999;
318  AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
319  if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
320  continue;
321  TrackMatch++; //this cluster has been matched, let's add it to the list
322  cluster_id.push_back(iClus);
323  Int_t ismatch = -1;
324  //now we need to find its matched corrected cluster if any
325  for (Int_t ic = 0;ic < nCluster2;ic++){
326  //we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
327  AliVCluster *emcluster2 = static_cast<AliVCluster*>(fCaloClusters2->At(ic));
328  if (! emcluster2)
329  continue;
330  if (! AcceptCluster(emcluster2))
331  continue;
332  if (! IsJetCluster(jet,ic,false))
333  continue;
334  TLorentzVector nPart;
335  cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
336  TLorentzVector nPart2;
337  emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
338  float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
339  if (R < 0.001){//this cluster stayed in the jet!
340  cluster_id2.push_back(ic);
341  ismatch = ic;
342  }
343  } // end of cluster2 loop
344  //we only get here if the track was matched to a cluster that hadn't been looked at before
345  if (ismatch < 0) //this cluster was entirely deleted
346  Esub+=cluster->E();
347  else{ //get the corrected cluster
348  AliVCluster *emcluster2temp = static_cast<AliVCluster*>(fCaloClusters2->At(ismatch));
349  Esub+=(cluster->E() - emcluster2temp->E());
350  }
351 
352  } // end of track loop
353  fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
354  fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
355  fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
356  fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
357  fHistNconvsPt[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
358  fHistNtvsPt[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
359  fHistNcvsPt[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
360  if (jet->MaxTrackPt()<5.0)
361  continue;
362  fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
363  fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
364  fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
365  fHistNconvsPtBias[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
366  fHistNtvsPtBias[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
367  fHistNcvsPtBias[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
368  fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
369  if (centbin == 0)
370  fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,jet->GetNumberOfTracks());
371  }
372  fHistNjetvsCent->Fill(fCent,NjetAcc);
373  return kTRUE;
374 }
375 
376 
377 
378 
379 
Double_t Area() const
Definition: AliEmcalJet.h:123
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
Bool_t AcceptTrack(AliVParticle *track, Int_t c=0) const
Bool_t fLocalInitialized
whether or not the task has been already initialized
Bool_t AcceptCluster(AliVCluster *clus, Int_t c=0) const
Cluster selection.
AliVTrack * GetTrack() const
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Bool_t IsJetTrack(AliEmcalJet *jet, Int_t itrack, Bool_t sorted=kFALSE) const
Bool_t IsJetCluster(AliEmcalJet *jet, Int_t iclus, Bool_t sorted=kFALSE) const
TClonesArray * fCaloClusters
!clusters
AliRhoParameter * GetRhoFromEvent(const char *name)
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:131
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
AliRhoParameter * fRho
! event rho
void ExecOnce()
Perform steps needed to initialize the analysis.
BeamType GetBeamType() const
Get beam type.
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:148
Double_t fCent
!event centrality
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:102
AliEmcalList * fOutput
!output list
AliVCluster * GetCluster() const
TClonesArray * fTracks
!tracks
virtual Int_t GetCentBin(Double_t cent) const
Double_t fVertex[3]
!event vertex
ClassImp(AliAnalysisTaskEmcalJetHadCorQA) AliAnalysisTaskEmcalJetHadCorQA
Bool_t fCreateHisto
whether or not create histograms
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
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.
Bool_t IsEMCAL() const
void UserCreateOutputObjects()
Main initialization function on the worker.
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
bool Bool_t
Definition: External.C:53
Double_t NEF() const
Definition: AliEmcalJet.h:141
Double_t fRhoVal
! event rho value, same for local rho
TH2F * fHistNEFvsPt[3]
number of jets versus Centrality
Float_t RelativePhi(Double_t mphi, Double_t vphi) const