1 #include "AliAnalysisTaskEmcalJetHadCorQA.h"
5 #include <TClonesArray.h>
11 #include <TLorentzVector.h>
12 #include <TParameter.h>
13 #include <TParticle.h>
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"
32 ClassImp(AliAnalysisTaskEmcalJetHadCorQA)
34 //________________________________________________________________________
35 AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() :
36 AliAnalysisTaskEmcalJet("jethadcor",kFALSE),
43 fHistNTMatchvsPtvsNtack0(0)
45 for (int i = 0;i<3;i++){
47 fHistNTMatchvsPt[i] = 0;
48 fHistNCMatchvsPt[i] = 0;
49 fHistHadCorvsPt[i] = 0;
50 fHistNEFvsPtBias[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;
61 // Default constructor.
63 SetMakeGeneralHistograms(kTRUE);
66 //________________________________________________________________________
67 AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
68 AliAnalysisTaskEmcalJet(name,kTRUE),
75 fHistNTMatchvsPtvsNtack0(0)
77 for (int i = 0;i<3;i++){
79 fHistNTMatchvsPt[i] = 0;
80 fHistNCMatchvsPt[i] = 0;
81 fHistHadCorvsPt[i] = 0;
82 fHistNEFvsPtBias[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;
93 SetMakeGeneralHistograms(kTRUE);
96 //________________________________________________________________________
97 void AliAnalysisTaskEmcalJetHadCorQA::UserCreateOutputObjects()
101 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
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);
105 for (int i = 0;i<3;i++){
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]);
151 fHistNTMatchvsPtvsNtack0 = new TH3F("NTMmatchvsPtvsNtrack0", "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);
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]);
160 fOutput->Add(fHistRhovsCent);
161 fOutput->Add(fHistNjetvsCent);
162 fOutput->Add(fHistNTMatchvsPtvsNtack0);
163 PostData(1, fOutput);
166 //________________________________________________________________________
168 Int_t AliAnalysisTaskEmcalJetHadCorQA::GetCentBin(Double_t cent) const
170 // Get centrality bin.
173 if (cent>=0 && cent<10)
175 else if (cent>=10 && cent<30)
177 else if (cent>=30 && cent<50)
184 //________________________________________________________________________
186 Float_t AliAnalysisTaskEmcalJetHadCorQA:: RelativePhi(Double_t mphi,Double_t vphi) const
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());
196 return dphi;//dphi in [-Pi, Pi]
199 //________________________________________________________________________
200 void AliAnalysisTaskEmcalJetHadCorQA::ExecOnce(){
201 AliAnalysisTaskEmcalJet::ExecOnce();
203 // AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
204 // if (fCaloName == "CaloClusters")
205 // am->LoadBranch("CaloClusters");
207 if (!fCalo2Name.IsNull() && !fCaloClusters2){
208 fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
209 if (!fCaloClusters2){
210 AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
211 fInitialized = kFALSE;
216 if (!fMCParticlesName.IsNull() && !fMCParticles){
217 fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
219 AliError(Form("%s: Could not retrieve MC Particles %s!",GetName(),fMCParticlesName.Data()));
220 fInitialized = kFALSE;
224 //TString fMCParticlesName;
225 // TClonesArray *fMCParticles;
227 // else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
228 // AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
229 // fCaloClusters2 = 0;
230 // fInitialized = kFALSE;
236 //________________________________________________________________________
237 Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
239 Int_t centbin = GetCentBin(fCent);
240 //for pp analyses we will just use the first centrality bin
241 if (GetBeamType()==0)
251 const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
252 const Int_t nTrack = fTracks->GetEntriesFast();
253 // const Int_t nMC = fMCParticles->GetEntriesFast();
255 TString fRhoScaledName = fRhoName;
256 fRho = GetRhoFromEvent(fRhoScaledName);
257 fRhoVal = fRho->GetVal();
258 if (GetBeamType()==0)
260 fHistRhovsCent->Fill(fCent,fRhoVal);
261 const Int_t Njets = fJets->GetEntriesFast();
266 for (Int_t iJets = 0; iJets < Njets; ++iJets) {
267 Int_t TrackMatch = 0;
268 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
275 if (jet->MaxTrackPt()>100)
277 if (! AcceptJet(jet))
280 vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
281 vector<Int_t> cluster_id2;
282 Double_t Esub = 0; //total E subtracted from the jet
283 Double_t jetPt = -500;
284 jetPt = jet->Pt()-jet->Area()*fRhoVal;
285 for (int i = 0;i<nTrack;i++){
286 AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(i));
289 if (!emctrack->IsEMCAL())
291 AliVTrack *track = (AliVTrack*)emctrack->GetTrack();
294 if (! AcceptTrack(track))
296 if (! IsJetTrack(jet,i,false))
298 Int_t iClus = track->GetEMCALcluster();
301 //we have the id of the matched cluster of a track constituent
302 // cout<<"track label for matched,accepted jet track is "<<track->GetLabel()<<endl;
303 bool ischecked = false;
304 for (UInt_t icid = 0;icid<cluster_id.size();icid++)
305 if (cluster_id[icid] == iClus)
306 ischecked = true; // we've already looked at this uncorrected cluster
308 continue; //no need to go further
309 AliEmcalParticle *emcluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
310 AliVCluster *cluster = emcluster->GetCluster();
313 if (! AcceptCluster(cluster))
315 Double_t etadiff = 999;
316 Double_t phidiff = 999;
317 AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
318 if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
320 TrackMatch++; //this cluster has been matched, let's add it to the list
321 cluster_id.push_back(iClus);
323 //now we need to find its matched corrected cluster if any
324 for (Int_t ic = 0;ic < nCluster2;ic++){
325 //we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
326 AliVCluster *emcluster2 = static_cast<AliVCluster*>(fCaloClusters2->At(ic));
329 if (! AcceptCluster(emcluster2))
331 if (! IsJetCluster(jet,ic,false))
333 TLorentzVector nPart;
334 cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
335 TLorentzVector nPart2;
336 emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
337 float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
338 if (R < 0.001){//this cluster stayed in the jet!
339 cluster_id2.push_back(ic);
342 } // end of cluster2 loop
343 //we only get here if the track was matched to a cluster that hadn't been looked at before
344 if (ismatch < 0) //this cluster was entirely deleted
346 else{ //get the corrected cluster
347 AliVCluster *emcluster2temp = static_cast<AliVCluster*>(fCaloClusters2->At(ismatch));
348 Esub+=(cluster->E() - emcluster2temp->E());
351 } // end of track loop
352 fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
353 fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
354 fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
355 fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
356 fHistNconvsPt[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
357 fHistNtvsPt[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
358 fHistNcvsPt[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
359 if (jet->MaxTrackPt()<5.0)
361 fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
362 fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
363 fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
364 fHistNconvsPtBias[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
365 fHistNtvsPtBias[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
366 fHistNcvsPtBias[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
367 fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
369 fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,jet->GetNumberOfTracks());
371 fHistNjetvsCent->Fill(fCent,NjetAcc);