]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Compilation with Root6
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetHF.cxx
CommitLineData
5e5a2b89 1//
2// Author: A Castro
3
4#include "AliAnalysisTaskEmcalJetHF.h"
5
6// general ROOT includes
7#include <TCanvas.h>
8#include <TChain.h>
9#include <TClonesArray.h>
10#include <TH1F.h>
11#include <TH2F.h>
12#include <TH3F.h>
13#include <THnSparse.h>
14#include <TList.h>
15#include <TLorentzVector.h>
16#include <TParameter.h>
17#include <TParticle.h>
18#include <TTree.h>
19#include <TVector3.h>
20#include <TObjArray.h>
21
22// AliROOT includes
23#include "AliAODEvent.h"
24#include "AliESDEvent.h"
25#include "AliAnalysisManager.h"
26#include "AliAnalysisTask.h"
27#include "AliCentrality.h"
28#include "AliEmcalJet.h"
29#include "AliAODJet.h"
30#include "AliVCluster.h"
31#include "AliVTrack.h"
32#include <AliVEvent.h>
33#include <AliVParticle.h>
34#include "AliRhoParameter.h"
35#include "AliEmcalParticle.h"
36
37// event handler (and pico's) includes
38#include <AliInputEventHandler.h>
39#include <AliVEventHandler.h>
40#include "AliESDInputHandler.h"
41#include "AliPicoTrack.h"
42#include "AliEventPoolManager.h"
43
44// PID includes
45#include "AliPIDResponse.h"
46#include "AliTPCPIDResponse.h"
47#include "AliESDpid.h"
48
49#include <AliInputEventHandler.h>
50#include <AliVEventHandler.h>
51
52// magnetic field includes
53#include "TGeoGlobalMagField.h"
54#include "AliMagF.h"
55
e742b9b0 56using std::cout;
57using std::endl;
58
5e5a2b89 59ClassImp(AliAnalysisTaskEmcalJetHF)
60
61//________________________________________________________________________
62AliAnalysisTaskEmcalJetHF::AliAnalysisTaskEmcalJetHF() :
63 AliAnalysisTaskEmcalJet("heavyF",kFALSE),
64 event(0),
65 fillHist(0),
08f52641 66 //isESD(0),
47d82e4b 67 doGlobalPID(0),
5e5a2b89 68 fCuts(0),
47d82e4b 69 fPhimin(-10), fPhimax(10),
70 fEtamin(-0.9), fEtamax(0.9),
71 fAreacut(0.0),
5e5a2b89 72 fJetHIpt(50.0),
73 fTrackPtCut(2.0),
74 fTrackEta(0.9),
08f52641 75 fTrkQAcut(0),
5e5a2b89 76 fPIDResponse(0x0), fTPCResponse(),
77 fEsdtrackCutsITSTPC(),
78 fEsdtrackCutsTPC(),
79 fEsdtrackCutsITS(),
80 fESD(0), fAOD(0),
81 fHistRhovsCent(0),
82 fHistJetPhi(0),
83 fHistCorJetPt(0), fHistJetPt(0),
84 fHistdEdx(0),
85 fHistdEdxvPt(0),
86 fHistClusE(0),
87 fHistEovPTracks(0),
88 fHistEovPvsPtTracks(0),
89 fHistPID(0), fHistPIDtpc(0), fHistPIDits(0), fHistPIDtof(0),
90 fHistnsigelectron(0),
91 fHistnSigElecPt(0),
92 fHistTrackPhivEta(0),
93 fHistClusterPhivEta(0),
94 fHistnJetTrackvnJetClusters(0),
47d82e4b 95 fhnPIDHF(0x0), fhnQA(0x0), fhnJetQA(0x0), fhnClusQA(0x0), fhnTrackQA(0x0), fhnGlobalPID(0x0)
5e5a2b89 96{
97 // Default constructor.
98 for (Int_t i = 0;i<6;++i){
99 fHistJetPtvsTrackPt[i] = 0;
100 fHistTrackPt[i] = 0;
101 fHistEP0[i] = 0;
102 fHistEP0A[i] = 0;
103 fHistEP0C[i] = 0;
104 fHistEPAvsC[i] = 0;
105 }
106
107 SetMakeGeneralHistograms(kTRUE);
108
109}
110
111//________________________________________________________________________
112AliAnalysisTaskEmcalJetHF::AliAnalysisTaskEmcalJetHF(const char *name) :
113 AliAnalysisTaskEmcalJet(name,kTRUE),
114 event(0),
115 fillHist(0),
08f52641 116 //isESD(0),
47d82e4b 117 doGlobalPID(0),
5e5a2b89 118 fCuts(0),
47d82e4b 119 fPhimin(-10), fPhimax(10),
120 fEtamin(-0.9), fEtamax(0.9),
121 fAreacut(0.0),
5e5a2b89 122 fJetHIpt(50.0),
123 fTrackPtCut(2.0),
124 fTrackEta(0.9),
08f52641 125 fTrkQAcut(0),
5e5a2b89 126 fPIDResponse(0x0), fTPCResponse(),
127 fEsdtrackCutsITSTPC(),
128 fEsdtrackCutsTPC(),
129 fEsdtrackCutsITS(),
130 fESD(0), fAOD(0),
131 fHistRhovsCent(0),
132 fHistJetPhi(0),
133 fHistCorJetPt(0), fHistJetPt(0),
134 fHistdEdx(0),
135 fHistdEdxvPt(0),
136 fHistClusE(0),
137 fHistEovPTracks(0),
138 fHistEovPvsPtTracks(0),
139 fHistPID(0), fHistPIDtpc(0), fHistPIDits(0), fHistPIDtof(0),
140 fHistnsigelectron(0),
141 fHistnSigElecPt(0),
142 fHistTrackPhivEta(0),
143 fHistClusterPhivEta(0),
144 fHistnJetTrackvnJetClusters(0),
47d82e4b 145 fhnPIDHF(0x0), fhnQA(0x0), fhnJetQA(0x0), fhnClusQA(0x0), fhnTrackQA(0x0), fhnGlobalPID(0x0)
5e5a2b89 146{
147 for (Int_t i = 0;i<6;++i){
148 fHistJetPtvsTrackPt[i] = 0;
149 fHistTrackPt[i] = 0;
150 fHistEP0[i] = 0;
151 fHistEP0A[i] = 0;
152 fHistEP0C[i] = 0;
153 fHistEPAvsC[i] = 0;
154 }
155 SetMakeGeneralHistograms(kTRUE);
156
157 DefineInput(0,TChain::Class());
158 DefineOutput(1, TList::Class());
159}
160
161//_______________________________________________________________________
162AliAnalysisTaskEmcalJetHF::~AliAnalysisTaskEmcalJetHF()
163{
164 // destructor
165 //
166 if (fOutput) {
167 delete fOutput;
168 fOutput = 0;
169 }
170}
171
172//________________________________________________________________________
173void AliAnalysisTaskEmcalJetHF::UserCreateOutputObjects()
174{
175 if (! fCreateHisto)
176 return;
177 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
178
179 fOutput = new TList();
180 fOutput->SetOwner(kTRUE);
181
182 fHistJetPhi = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2*TMath::Pi(),2*TMath::Pi());
183 fHistJetPt = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
184 fOutput->Add(fHistJetPhi);
185 fOutput->Add(fHistJetPt);
186
187 fillHist = 1;
188
189 fHistClusE = new TH1F("NumberClustersvsEnergy","NumberClustersvsEnergy", 500, 0, 10);
190
191
192 if(fillHist>0){
193 fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 400, 0, 400);
194 fHistCorJetPt = new TH1F("NjetvsCorrJetPt", "NjetvsCorrJetPt", 300, -100, 200);
195 fHistdEdx = new TH1F("dEdxSignal", "dEdxSignal", 500, 0, 500);
196 fHistdEdxvPt = new TH2F("dEdxvPt", "dEdxvPt", 200, 0, 100,500, 0 ,500);
197 fHistEovPTracks = new TH1F("EovPTracks","EovPTracks",200,0.0,2.0);
198 fHistEovPvsPtTracks = new TH2F("E/p_vs_Pt","E/p_vs_Pt",200, 0 ,100 ,50, 0, 2);
199 fHistPID = new TH1F("fHistPID", "Detector PID", 8, 0, 8);
200 fHistPIDtpc = new TH1F("fHistPIDtpc", "TPC pid", 4, 0, 4);
201 fHistPIDits = new TH1F("fHistPIDits", "ITS pid", 4, 0, 4);
202 fHistPIDtof = new TH1F("fHistPIDtof", "TOF pid", 4, 0, 4);
203 fHistnsigelectron = new TH1F("nsig_elec_TPC","nsig_elec_TPC",500,-10,10);
204 fHistnSigElecPt = new TH2F("nsig_v_pt(TPC)","nsig_v_pt(TPC)",200,0,100,100,-10,10);
205 fHistTrackPhivEta = new TH2F("TrackPhi_v_Eta","TrackPhiEta",64,-1.6,1.6,72,1,2*TMath::Pi());
206 fHistClusterPhivEta = new TH2F("ClusterPhi_v_Eta","ClusterPhiEta",64,-1.6,1.6,72,1,2*TMath::Pi());
207 fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
208
209 // PT bins used to be (2000, -100, 300)
210 TString name;
211 TString title;
212
213 // creating centrality dependent histos that don't involve Global Rho
214 for (Int_t i = 0;i<6;++i){
215 name = TString(Form("JetPtvsTrackPt_%i",i));
216 title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
217 fHistJetPtvsTrackPt[i] = new TH2F(name,title, 500, -100, 400, 100,0,100);
218 fOutput->Add(fHistJetPtvsTrackPt[i]);
219
220 name = TString(Form("TrackPt_%i",i));
221 title = TString(Form("Track pT cent bin %i",i));
222 fHistTrackPt[i] = new TH1F(name,title,400,0,200);
223 fOutput->Add(fHistTrackPt[i]);
224
225 name = TString(Form("EP0_%i",i));
226 title = TString(Form("EP VZero cent bin %i",i));
227 fHistEP0[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
228 fOutput->Add(fHistEP0[i]);
229
230 name = TString(Form("EP0A_%i",i));
231 title = TString(Form("EP VZero cent bin %i",i));
232 fHistEP0A[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
233 fOutput->Add(fHistEP0A[i]);
234
235 name = TString(Form("EP0C_%i",i));
236 title = TString(Form("EP VZero cent bin %i",i));
237 fHistEP0C[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
238 fOutput->Add(fHistEP0C[i]);
239
240 name = TString(Form("EPAvsC_%i",i));
241 title = TString(Form("EP VZero cent bin %i",i));
242 fHistEPAvsC[i] = new TH2F(name,title,144,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
243 fOutput->Add(fHistEPAvsC[i]);
244
245 }
246
247 fOutput->Add(fHistRhovsCent);
248 fOutput->Add(fHistCorJetPt);
249 fOutput->Add(fHistdEdx);
250 fOutput->Add(fHistdEdxvPt);
251 fOutput->Add(fHistEovPTracks);
252 fOutput->Add(fHistEovPvsPtTracks);
253 fOutput->Add(fHistPID);
254 fOutput->Add(fHistPIDtpc);
255 fOutput->Add(fHistPIDits);
256 fOutput->Add(fHistPIDtof);
257 fOutput->Add(fHistnsigelectron);
258 fOutput->Add(fHistnSigElecPt);
259 fOutput->Add(fHistTrackPhivEta);
260 fOutput->Add(fHistClusterPhivEta);
261 fOutput->Add(fHistnJetTrackvnJetClusters);
262
263 }//Fill Histograms
264
265 fOutput->Add(fHistClusE);
266
267
268 // ****************************** PID *****************************************************
269 // set up PID handler
270 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
271 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
272 if(!inputHandler) {
273 AliFatal("Input handler needed");
274 }
275
276 // PID response object
277 //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
278 // inputHandler->CreatePIDResponse(fIsMC); // needed to create object, why though?
279 fPIDResponse = inputHandler->GetPIDResponse();
280 if (!fPIDResponse) {
281 AliError("PIDResponse object was not created");
282 }
08f52641 283 // ****************************************************************************************
284 UInt_t bitcoded = 0; // bit coded, see GetDimParamsPID() below
285 bitcoded = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
5e5a2b89 286 fhnPIDHF = NewTHnSparseDHF("fhnPIDHF", bitcoded);
08f52641 287
288 UInt_t bitcoded1 = 0; // bit coded, see GetDimParamsPID() below
289 bitcoded1 = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4;
290 fhnJetQA = NewTHnSparseDJetQA("fhnJetQA", bitcoded1);
291
292 UInt_t bitcoded2 = 0; // bit coded, see GetDimParamsPID() below
293 bitcoded2 = 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15;
294 fhnClusQA = NewTHnSparseDJetQA("fhnClusQA", bitcoded2);
295
296 UInt_t bitcoded3 = 0; // bit coded, see GetDimParamsPID() below
297 bitcoded3 = 1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14;
298 fhnTrackQA = NewTHnSparseDJetQA("fhnTrackQA", bitcoded3);
299
5e5a2b89 300 cout << "_______________Created Sparse__________________" << endl;
08f52641 301
5e5a2b89 302 fOutput->Add(fhnPIDHF);
08f52641 303 fOutput->Add(fhnJetQA);
304 fOutput->Add(fhnClusQA);
305 fOutput->Add(fhnTrackQA);
47d82e4b 306 fOutput->Add(fhnGlobalPID);
08f52641 307
5e5a2b89 308
309 PostData(1, fOutput);
310}
311
312//________________________________________________________
313void AliAnalysisTaskEmcalJetHF::ExecOnce()
314{
315 // Initialize the analysis
316 AliAnalysisTaskEmcalJet::ExecOnce();
317
318} // end of ExecOnce
319
320//________________________________________________________________________
321Bool_t AliAnalysisTaskEmcalJetHF::Run()
322{
323 // check to see if we have any tracks
324 if (!fTracks) return kTRUE;
325 if (!fJets) return kTRUE;
08f52641 326
08f52641 327 // what kind of event do we have: AOD or ESD?
47d82e4b 328 Bool_t useAOD;
08f52641 329 if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
330 else useAOD = kFALSE;
331
332 // if we have ESD event, set up ESD object
333 if(!useAOD){
334 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
335 if (!fESD) {
336 AliError(Form("ERROR: fESD not available\n"));
337 return kTRUE;
5e5a2b89 338 }
5e5a2b89 339 }
340
341 // if we have AOD event, set up AOD object
342 if(useAOD){
343 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
344 if(!fAOD) {
345 AliError(Form("ERROR: fAOD not available\n"));
346 return kTRUE;
347 }
348 }
5e5a2b89 349
350 // get magnetic field info for DCA
351 Double_t MagF = fESD->GetMagneticField();
352 Double_t MagSign = 1.0;
353 if(MagF<0)MagSign = -1.0;
354 // set magnetic field
355 if (!TGeoGlobalMagField::Instance()->GetField()) {
356 AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
357 TGeoGlobalMagField::Instance()->SetField(field);
358 }
359
360 // get centrality bin
361 Int_t centbin = GetCentBin(fCent);
362 //for pp analyses we will just use the first centrality bin
363 if (centbin == -1) centbin = 0;
364
365 // get vertex information
366 Double_t fvertex[3]={0,0,0};
367 InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
368 //Double_t zVtx=fvertex[2];
369
370 // create pointer to list of input event
371 TList *list = InputEvent()->GetList();
372 if(!list) {
373 AliError(Form("ERROR: list not attached\n"));
374 return kTRUE;
375 }
376
377 // background density
378 fRhoVal = fRho->GetVal();
379
380 // initialize TClonesArray pointers to jets and tracks
381 TClonesArray *jets = 0;
382 TClonesArray *tracks = 0;
383 TClonesArray *clusters = 0;
384
385 // get Jets object
386 jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
387 if(!jets){
388 AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
389 return kTRUE;
390 } // verify existence of jets
391
392 // get Tracks object
393 tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
394 if (!tracks) {
395 AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
396 return kTRUE;
397 } // verify existence of tracks (fTracksName.Data())
398
399 //get Clusters object
400 clusters = dynamic_cast<TClonesArray*>(list->FindObject(fCaloClusters));
401 if (!clusters){
402 AliError(Form("Pointer to tracks %s == 0",fCaloClusters->GetName()));
403 return kTRUE;
404 } //verify cluster existence
08f52641 405
406
407 //const Int_t nclus = clusters->GetEntries();
408 /* Globalcluster
409 for(int icluster = 0; icluster < nclus; icluster++) {
410 AliVCluster* andyCluster = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
411 if (!icluster) continue;
412 }
413 */
5e5a2b89 414 // get all emcal clusters
415 TRefArray* caloClusters = new TRefArray();
416 fESD->GetEMCALClusters( caloClusters );
417
418 //TObjArray* listcuts = fEsdtrackCutsTPC->GetAcceptedTracks(fESD);
419 //Int_t nGoodTracks = list->GetEntries();
420 Int_t nCluster = caloClusters->GetEntries();
421
422 // event plane histograms filled
423 if(fillHist>0) fHistEP0[centbin]->Fill(fEPV0);
424 if(fillHist>0) fHistEP0A[centbin]->Fill(fEPV0A);
425 if(fillHist>0) fHistEP0C[centbin]->Fill(fEPV0C);
426 if(fillHist>0) fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
47d82e4b 427
428 event++;
429 cout<<"Event #: "<<event<<endl;
430
5e5a2b89 431 // get number of jets and tracks
432 const Int_t Njets = jets->GetEntries();
433 const Int_t Ntracks = tracks->GetEntries();
434 //const Int_t Nclusters = clusters->GetEntries();
435 if(Ntracks<1) return kTRUE;
436 if(Njets<1) return kTRUE;
437 //if(Nclusters<1) return kTRUE;
47d82e4b 438
439 if(doGlobalPID){
440 // loop over Globap tracks in the event for PID: Runs ON ESD's
5e5a2b89 441 for (int i = 0;i<Ntracks;i++){
442 //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
47d82e4b 443 AliVTrack *trackGlobal = static_cast<AliVTrack*>(fTracks->At(i));
08f52641 444 //AliESDtrack* track = fESD->GetTrack(i);
47d82e4b 445 if (! trackGlobal) {
5e5a2b89 446 AliError(Form("Couldn't get ESD track %d\n", i));
447 continue;
448 }
449
450
451 // apply track cuts
47d82e4b 452 if(TMath::Abs(trackGlobal->Eta())>fTrackEta) continue;
453 if (trackGlobal->Pt()<0.15) continue;
5e5a2b89 454
455 // initialize track info
47d82e4b 456 Double_t dEdxGlobal = -99;
457 Double_t trackphiGlobal = -99;
458 Double_t trackptGlobal = 0;
459 Double_t p = trackGlobal->P();
460 Double_t fClsEGlobal = -99;
461 Double_t EovPGlobal = -99;
5e5a2b89 462
463 // distance of closest approach
47d82e4b 464 Float_t DCAxyGlobal = -999;
465 Float_t DCAzGlobal = -999;
466 //Double_t DCAXYGlobal = -999;
467 //Double_t DCAZGlobal = -999;
5e5a2b89 468
469 // track info
47d82e4b 470 trackGlobal->GetTPCsignal();
471 //trackGlobal->GetImpactParameters(DCAxyGlobal, DCAzGlobal);
472 dEdxGlobal = trackGlobal->GetTPCsignal();
473 trackptGlobal = trackGlobal->Pt();
474 trackphiGlobal = trackGlobal->Phi();
5e5a2b89 475
476 // fill track histo's
47d82e4b 477 //if(fillHist>0) fHistdEdx->Fill(dEdx);
478 //if(fillHist>0) fHistdEdxvPt->Fill(trackpt,dEdx);
479 //if(fillHist>0) fHistTrackPt[centbin]->Fill(track->Pt());
5e5a2b89 480
481 // clusters
47d82e4b 482 Int_t nMatchClusGlobal = -1;
483 AliESDCaloCluster *matchedClusGlobal = NULL;
5e5a2b89 484
485 //////////////////////////////////////////////////////////////////////////
08f52641 486
5e5a2b89 487 // cut on 1.5 GeV for EMCal Cluster //if(track->GetEMCALcluster()<0 || pt<1.5) continue;
488 //////////////////////////////////////////////////////////////////////////
489 // particles in TOF
47d82e4b 490 //Double_t nSigmaPion_TOFGlobal, nSigmaProton_TOFGlobal, nSigmaKaon_TOFGlobal = -1.;
5e5a2b89 491
492 // misc quantities
47d82e4b 493 //Double_t TOFsigGlobal = -1.;
494 //Double_t nClustersTPCGlobal = -1;
495 //Int_t chargeGlobal = 0;
496 //Int_t trackCutsGlobal = 0; // set to 0 to get to run
497 //Int_t myPIDGlboal = 0; // set to 0 because myPID is to compare with MC unless I hard code in cuts
498 Double_t nSigmaElectron_TPCGlobal = -999;
499 Double_t nSigmaElectron_TOFGlobal = -999;
500 //Double_t nSigmaElectron_EMCALGlobal = -999;
5e5a2b89 501
502 // get clusters
47d82e4b 503 //Int_t clusiGlobal = trackGlobal->GetEMCALcluster();
504 //nClustersTPCGlobal = trackGlobal->GetTPCclusters(0);
5e5a2b89 505
506 AliESDtrack *trackESD = fESD->GetTrack(Ntracks);
507
47d82e4b 508 nSigmaElectron_TPCGlobal = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
509 nSigmaElectron_TOFGlobal = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron);
5e5a2b89 510
511 //for EMCAL
47d82e4b 512 nMatchClusGlobal = trackGlobal->GetEMCALcluster();
513 if(nMatchClusGlobal > 0){
514 matchedClusGlobal = (AliESDCaloCluster*)fESD->GetCaloCluster(nMatchClusGlobal);
5e5a2b89 515 //AliESDCaloCluster* matchedClus = fESD->GetCaloCluster(clusi);
47d82e4b 516 //double eop2Global = -1;
517 //double ssGl[4]={0.,0.,0.,0.};
5e5a2b89 518 //Double_t nSigmaEop = fPID->GetPIDResponse()-m >NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
519
47d82e4b 520 fClsEGlobal = matchedClusGlobal->E();
521 EovPGlobal = fClsEGlobal/p;
522 //nSigmaElectron_EMCALGlobal = fPIDResponse->NumberOfSigmasEMCAL(trackESD,AliPID::kElectron,eop2,ss);
523 //if(fillHist>0) fHistEovPTracks->Fill(EovP);
524 //if(fillHist>0) fHistEovPvsPtTracks->Fill(trackpt,EovP);
5e5a2b89 525 }
526
47d82e4b 527 //if(fillHist>0) fHistnsigelectron->Fill(nSigmaElectron_TPC);
528 //if(fillHist>0) fHistnSigElecPt->Fill(trackpt,nSigmaElectron_TPC);
5e5a2b89 529
530 // extra attempt
47d82e4b 531 AliVEvent *eventQA=InputEvent();
532 if (!eventQA||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
5e5a2b89 533 //////////////////////////////////////////////////////////////////////////
08f52641 534 // cut on 1.5 GeV for EMCal Cluster
47d82e4b 535 if(trackGlobal->Pt() < fTrackPtCut) continue;
5e5a2b89 536
47d82e4b 537 Double_t HF_tracks[12] = {fCent, trackGlobal->Pt(), trackGlobal->P() ,trackGlobal->Eta(), trackGlobal->Phi(), EovPGlobal, DCAxyGlobal, DCAzGlobal, dEdxGlobal, nSigmaElectron_TPCGlobal, nSigmaElectron_TOFGlobal, 0};//nSigmaElectron_EMCAL};
538 fhnGlobalPID->Fill(HF_tracks); // fill Sparse Histo with trigger entries
5e5a2b89 539
540 } // Loop over tracks
47d82e4b 541 }//PID Switch
08f52641 542
47d82e4b 543 // Start Jet Analysis
5e5a2b89 544 // initialize jet parameters
545 Int_t ijethi=-1;
546 Double_t highestjetpt=0.0;
547
08f52641 548 // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
5e5a2b89 549 for (Int_t ijet = 0; ijet < Njets; ijet++){
550 // get our jets
551 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
552 if (!jet) continue;
08f52641 553
47d82e4b 554 // apply jet cuts
555 if(!AcceptMyJet(jet)) continue;
556
e742b9b0 557 Double_t JetQA[5] = {0, static_cast<Double_t>(jet->GetNumberOfTracks()), static_cast<Double_t>(jet->GetNumberOfClusters()),jet->Eta(), jet->Phi()};
08f52641 558 fhnJetQA->Fill(JetQA);
559
560
561 // Loop over clusters for JetQA
562 for(int iCluster = 0; iCluster < nCluster; iCluster++) {
563 AliVCluster* caloCluster = static_cast<AliVCluster*>(fCaloClusters->At(iCluster));
564 //AliVCluster* caloCluster = (AliVCluster* )caloClusters->At(jet->GetNumberOfClusters());
565 //AliESDCaloCluster* caloCluster = (AliESDCaloCluster* )caloClusters->At(iCluster);
566 //AliESDCaloCluster* clus = fESD->GetCaloCluster(iclus);
567 if (!caloCluster){
568 AliError(Form("ERROR: Could not get cluster %d", iCluster));
569 continue;
570 }
571 if(!IsJetCluster(jet,iCluster,kFALSE)) continue ;
572 Float_t pos[3];
573 caloCluster->GetPosition(pos); // Get cluster position
574 TVector3 cp(pos);
575 Double_t NtrMatched = -999.0;
576 NtrMatched = caloCluster->GetNTracksMatched();
577
578 //loop over tracks for Jet QA
579 for(int itrack = 0; itrack < NtrMatched; itrack++){
580 //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
581 AliVTrack *trackcluster = static_cast<AliVTrack*>(fTracks->At(itrack));
582 //AliESDtrack* track = fESD->GetTrack(i);
583 if (! trackcluster) {
584 AliError(Form("Couldn't get ESD track %d\n", itrack));
585 continue;
586 }
587 if(!IsJetTrack(jet,itrack,kFALSE)) continue;
e742b9b0 588 Double_t ClusQA[6] = {caloCluster->E(),cp.PseudoRapidity() ,cp.Phi(), static_cast<Double_t>(caloCluster->IsEMCAL()), NtrMatched, trackcluster->Phi()};
08f52641 589 fhnClusQA->Fill(ClusQA); //,1./Njets); // fill Sparse Histo with trigger entries
590 }//loop over tracks for JetQA
591
592 }//loop over clusters for JetQA
593
594 // loop over tracks in the event
595 for (int iTrack = 0; iTrack<jet->GetNumberOfTracks(); iTrack++){ //loop over tracks in jet for TrackQA
596 //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
597 AliVTrack *jetTrack = static_cast<AliVTrack*>(fTracks->At(iTrack));
598 //AliESDtrack* track = fESD->GetTrack(i);
599 if (! jetTrack) {
600 AliError(Form("Couldn't get ESD track %d\n", iTrack));
601 continue;
602 }
603 if(!IsJetTrack(jet,iTrack,kFALSE)) continue;
e742b9b0 604 Double_t trackQA[5] = {jetTrack->Pt(), jetTrack->Eta(), jetTrack->Phi(), static_cast<Double_t>(jetTrack->IsEMCAL()), static_cast<Double_t>(jetTrack->GetEMCALcluster())};
08f52641 605 fhnTrackQA->Fill(trackQA); //,1./Njets);
606
607 }//track loop for TrackQA
5e5a2b89 608
609 if(highestjetpt<jet->Pt()){
610 ijethi=ijet;
611 highestjetpt=jet->Pt();
612 }
613 } // end of looping over jets
614
615 // loop over jets in the event and make appropriate cuts
616 for (Int_t iJets = 0; iJets < Njets; ++iJets) {
617 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
618 if (!jet) // see if we have a jet
619 continue;
620
621 // phi of jet, constrained to 1.6 < Phi < 2.94
622 float jetphi = jet->Phi(); // phi of jet
623
47d82e4b 624 // apply jet cuts
625 if(!AcceptMyJet(jet)) continue;
626
5e5a2b89 627 Int_t JetClusters = jet->GetNumberOfClusters();
628 Int_t JetTracks = jet -> GetNumberOfTracks();
629 fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
08f52641 630
5e5a2b89 631 // Initializations and Calculations
632 //Double_t jeteta = jet->Eta(); // ETA of jet
633 Double_t jetptraw = jet->Pt(); // raw pT of jet
634 Double_t jetPt = -500; // initialize corr jet pt LOCAL
635 Double_t jetarea = -500; // initialize jet area
636 jetarea = jet->Area(); // jet area
637 jetPt = jet->Pt() - jetarea*fRhoVal; // semi-corrected pT of jet from GLOBAL rho value
638
639 // make histo's
640 if(fillHist>0) fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
641 fHistJetPhi->Fill(jetphi);
642 if(fillHist>0) fHistCorJetPt->Fill(jetPt);
643 fHistJetPt->Fill(jetptraw);
644
47d82e4b 645 if(jet->Pt() > fJetHIpt) {
5e5a2b89 646
5e5a2b89 647 //loop over clusters
648 //for (int i = 0; i < njetclusters; i++){
649 for(int iCluster = 0; iCluster < nCluster; iCluster++) {
08f52641 650 AliVCluster* caloCluster = static_cast<AliVCluster*>(fCaloClusters->At(iCluster));
651 //AliVCluster* caloCluster = (AliVCluster* )caloClusters->At(iCluster);
5e5a2b89 652 //AliESDCaloCluster* caloCluster = (AliESDCaloCluster* )caloClusters->At(iCluster);
653 //AliESDCaloCluster* clus = fESD->GetCaloCluster(iclus);
654 if (!caloCluster){
655 AliError(Form("ERROR: Could not get cluster %d", iCluster));
656 continue;
657 }
658
659 //if (!caloCluster -> IsEMCAL()) continue; //Check that Cluster is EMCal Cluster
660 if(! IsJetCluster(jet, iCluster, kFALSE)) continue;
661
08f52641 662 AliESDtrack *track = 0;
08f52641 663 if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
664 track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
5e5a2b89 665 Double_t fclusE = -999;
666 fclusE = caloCluster->E();
08f52641 667
5e5a2b89 668 if (fclusE<0.) continue; //Check that cluster has positive energy
669
670 fHistClusE->Fill(fclusE);
671 //Int_t cluslabel = caloCluster->GetID();
672 //Int_t nmatched = caloCluster->GetNTracksMatched();
673
674 Float_t pos[3];
675 caloCluster->GetPosition(pos); // Get cluster position
676 TVector3 cp(pos);
5e5a2b89 677 if(fillHist>0) fHistClusterPhivEta->Fill(cp.PseudoRapidity(),cp.Phi());
678
679 Int_t trackMatchedindex=caloCluster->GetTrackMatchedIndex();
5e5a2b89 680 if(trackMatchedindex<0)continue; // Make sure we don't have a bad index
08f52641 681
682 if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
683 track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
684 Double_t NtrMatched = -999;
685 NtrMatched = caloCluster->GetNTracksMatched();
686 for(int itrack = 0; itrack < NtrMatched; itrack++){ // Loop over tracks matched to clusters from jets
687 //AliVParticle *trackCluster = static_cast<AliVParticle*>(fTracks->At(i));
688 AliVTrack *trackCluster = static_cast<AliVTrack*>(fTracks->At(itrack));
689 //AliESDtrack* trackCluster = fESD->GetTrack(itrack);
690 if (! trackCluster) {
691 AliError(Form("Couldn't get ESD track %d\n", itrack));
692 continue;
693 }
694 if(!IsJetTrack(jet,itrack,kFALSE)) continue; // Check that track is still part of ith jet
5e5a2b89 695
08f52641 696 //if (trackCluster->Phi()>3.2 || trackCluster->Phi()<1.4)cout <<"Out of Range Track! Track Phi: " << trackCluster->Phi() << endl;
697
5e5a2b89 698 //if(track->Pt() < 4.0) continue;
5e5a2b89 699 // initialize track info
700 Double_t dEdx = -99;
701 //Double_t trackphi = -99;
702 Double_t trackpt = 0;
08f52641 703 Double_t p = trackCluster->P();
5e5a2b89 704 Double_t EovP = -99;
705
706 // distance of closest approach
707 Float_t DCAxy = -999;
708 Float_t DCAz = -999;
709 //Double_t DCAXY = -999;
710 //Double_t DCAZ = -999;
711
712 // track info
08f52641 713 //trackCluster->GetImpactParameters(DCAxy, DCAz);
714 trackpt = trackCluster->Pt();
5e5a2b89 715 //trackphi = track->Phi();
716
717 // fill track histo's
08f52641 718 if(fillHist>0) fHistTrackPhivEta->Fill(trackCluster->Eta(),trackCluster->Phi());
5e5a2b89 719 if(fillHist>0) fHistdEdx->Fill(dEdx);
720 if(fillHist>0) fHistdEdxvPt->Fill(trackpt,dEdx);
08f52641 721 if(fillHist>0) fHistTrackPt[centbin]->Fill(trackCluster->Pt());
722
5e5a2b89 723 Double_t nSigmaElectron_TPC = -999;
724 Double_t nSigmaElectron_TOF = -999;
725 Double_t nSigmaElectron_EMCAL = -999;
08f52641 726
47d82e4b 727 if(!useAOD){
5e5a2b89 728
729 AliESDtrack *trackESD = fESD->GetTrack(Ntracks);
730
08f52641 731 dEdx = trackESD->GetTPCsignal();
732 trackESD->GetImpactParameters(DCAxy, DCAz);
5e5a2b89 733 nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
734 nSigmaElectron_TOF = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron);
47d82e4b 735 }
736 if (useAOD) {
5e5a2b89 737 AliAODTrack *trackAOD = fAOD->GetTrack(Ntracks);
738
739 // get detector signals
740 dEdx = trackAOD->GetTPCsignal();
47d82e4b 741
742 //trackAOD->GetImpactParameters(DCAxy,DCAz);
5e5a2b89 743 nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(trackAOD,AliPID::kElectron);
744 nSigmaElectron_TOF = fPIDResponse->NumberOfSigmasTOF(trackAOD,AliPID::kElectron);
745
746 } // end of AOD pid
47d82e4b 747
5e5a2b89 748 //double eop2 = -1;
749 //double ss[4]={0.,0.,0.,0.};
750 EovP = fclusE/p;
751 //nSigmaElectron_EMCAL = fPIDResponse->NumberOfSigmasEMCAL(trackESD,AliPID::kElectron,eop2,ss);
752 if(fillHist>0) fHistEovPTracks->Fill(EovP);
753 if(fillHist>0) fHistEovPvsPtTracks->Fill(trackpt,EovP);
08f52641 754 Double_t HF_tracks[12] = {fCent, trackCluster->Pt(), trackCluster->P() ,trackCluster->Eta(), trackCluster->Phi(), EovP, 0/*DCAxy*/, 0/*DCAz*/, dEdx, nSigmaElectron_TPC, nSigmaElectron_TOF, nSigmaElectron_EMCAL};
5e5a2b89 755 fhnPIDHF->Fill(HF_tracks); // fill Sparse Histo with trigger entries
08f52641 756
5e5a2b89 757 //Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();//find the index of the matched track. This by default returns the best match
5e5a2b89 758 //AliESDtrack *track = event->GetTrack(trackMatchedIndex);
759 //if this is a good track, accept track will return true. The track matched is good, so not track matched is false
760 //Int_t matched = fESD->AcceptTrack(track);//If the track is bad, don't count it. By default even bad tracks are accepted
761
08f52641 762 } //loop over tracks matched to clusters in jet
763
5e5a2b89 764 } // loop over jet clusters
765
766
47d82e4b 767 } // highest pt jet cut
5e5a2b89 768 } // LOOP over JETS in event
769
5e5a2b89 770 return kTRUE;
771
772}
773//________________________________________________________________________
774void AliAnalysisTaskEmcalJetHF::Terminate(Option_t *)
775{
776 cout<<"###########################"<<endl;
777 cout<<"#### Task Finished ####"<<endl;
778 cout<<"###########################"<<endl;
779 cout<<"###########################"<<endl;
780} // end of terminate
781
47d82e4b 782//________________________________________________________________________
783Int_t AliAnalysisTaskEmcalJetHF::AcceptMyJet(AliEmcalJet *jet) {
784 //applies all jet cuts except pt
785 if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
786 if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
787 if (jet->Area()<fAreacut) return 0;
788 // prevents 0 area jets from sneaking by when area cut == 0
789 if (jet->Area()==0) return 0;
790 //exclude jets with extremely high pt tracks which are likely misreconstructed
791 if(jet->MaxTrackPt()>100) return 0;
792
793 //passed all above cuts
794 return 1;
795}
796
5e5a2b89 797//________________________________________________________________________
798Int_t AliAnalysisTaskEmcalJetHF::GetCentBin(Double_t cent) const
799{ // Get centrality bin.
800
801 Int_t centbin = -1;
802 if (cent>=0 && cent<10)
803 centbin = 0;
804 else if (cent>=10 && cent<20)
805 centbin = 1;
806 else if (cent>=20 && cent<30)
807 centbin = 2;
808 else if (cent>=30 && cent<40)
809 centbin = 3;
810 else if (cent>=40 && cent<50)
811 centbin = 4;
812 else if (cent>=50 && cent<90)
813 centbin = 5;
814 return centbin;
815}
816//____________________________________________________________________________________________
817THnSparse* AliAnalysisTaskEmcalJetHF::NewTHnSparseDHF(const char* name, UInt_t entries)
818{
819 // generate new THnSparseD PID, axes are defined in GetDimParams()
820 Int_t count = 0;
821 UInt_t tmp = entries;
822 while(tmp!=0){
823 count++;
824 tmp = tmp &~ -tmp; // clear lowest bit
825 }
826
827 TString hnTitle(name);
828 const Int_t dim = count;
829 Int_t nbins[dim];
830 Double_t xmin[dim];
831 Double_t xmax[dim];
832
833 Int_t i=0;
834 Int_t c=0;
835 while(c<dim && i<32){
836 if(entries&(1<<i)){
837
838 TString label("");
839 GetDimParamsHF(i, label, nbins[c], xmin[c], xmax[c]);
840 hnTitle += Form(";%s",label.Data());
841 c++;
842 }
843
844 i++;
845 }
846 hnTitle += ";";
847
848 return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
08f52641 849} // end of NewTHnSparseF PID
850
851THnSparse* AliAnalysisTaskEmcalJetHF::NewTHnSparseDJetQA(const char* name, UInt_t entries)
852{
853 // generate new THnSparseD JetQA, axes are defined in GetDimParamsJetQA()
854 Int_t count = 0;
855 UInt_t tmp = entries;
856 while(tmp!=0){
857 count++;
858 tmp = tmp &~ -tmp; // clear lowest bit
859 }
860
861 TString hnTitle(name);
862 const Int_t dim = count;
863 Int_t nbins[dim];
864 Double_t xmin[dim];
865 Double_t xmax[dim];
866
867 Int_t i=0;
868 Int_t c=0;
869 while(c<dim && i<32){
870 if(entries&(1<<i)){
871
872 TString label("");
873 GetDimParamsJetQA(i, label, nbins[c], xmin[c], xmax[c]);
874 hnTitle += Form(";%s",label.Data());
875 c++;
876 }
877
878 i++;
879 }
880 hnTitle += ";";
881
882 return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
883} // end of NewTHnSparseF JetQA
884
5e5a2b89 885
886void AliAnalysisTaskEmcalJetHF::GetDimParamsHF(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
887{
888 // stores label and binning of axis for THnSparse
889 const Double_t pi = TMath::Pi();
890
891 switch(iEntry){
892
893 case 0:
894 label = "V0 centrality (%)";
895 nbins = 10;
896 xmin = 0.;
897 xmax = 100.;
898 break;
899
900 case 1:
901 label = "Track p_{T}";
902 nbins = 750;
903 xmin = 0.;
904 xmax = 75.;
905 break;
906
907 case 2:
908 label = "Track p";
909 nbins = 750;
910 xmin = 0.;
911 xmax = 75.;
912 break;
913
914 case 3:
915 label = "Track Eta";
916 nbins = 64;
917 xmin = -1.6;
918 xmax = 1.6;
919 break;
920
921 case 4:
922 label = "Track Phi";
923 nbins = 72;
924 xmin = 0;
925 xmax = 2*pi;
926 break;
927
928 case 5:
929 label = "E/p of track";
930 nbins = 400;
931 xmin = 0;
932 xmax = 4.0;
933 break;
934
935 case 6:
936 label = "DCA xy";
937 nbins = 20;
938 xmin = -10;
939 xmax = 10;
940 break;
941
942 case 7:
943 label = "DCA z";
944 nbins = 20;
945 xmin = -10;
946 xmax = 10;
947 break;
948
949 case 8:
950 label = "dEdX of track - TPC";
951 nbins = 300;
952 xmin = 0;
953 xmax = 300;
954 break;
955
956 case 9:
957 label = "nSigma electron TPC";
958 nbins = 50;
959 xmin = -5;
960 xmax = 5;
961 break;
962
963 case 10:
964 label = "nSigma electron TOF";
965 nbins = 50;
966 xmin = -5;
967 xmax = 5;
968 break;
969
970 case 11:
971 label = "nSigma electron Emcal";
972 nbins = 50;
973 xmin = -5;
974 xmax = 5;
975 break;
976
977
978 } // end of switch
08f52641 979} // end of getting dim-params
980
981void AliAnalysisTaskEmcalJetHF::GetDimParamsJetQA(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
982{
983 // stores label and binning of axis for THnSparse
984 const Double_t pi = TMath::Pi();
985
986 switch(iEntry){
987
988 case 0:
989 label = "number of Jets in Event";
990 nbins = 2;
991 xmin = 0.;
992 xmax = 1.;
993 break;
994
995 case 1:
996 label = "number of Clusters in a Jet";
997 nbins = 100;
998 xmin = 0.;
999 xmax = 100.;
1000 break;
1001
1002 case 2:
1003 label = "number of Tracks in a Jet";
1004 nbins = 100;
1005 xmin = 0.;
1006 xmax = 100.;
1007 break;
1008
1009 case 3:
1010 label = "Jet Eta";
1011 nbins = 24;
1012 xmin = -1.2;
1013 xmax = 1.2;
1014 break;
1015
1016 case 4:
1017 label = "Jet Phi";
1018 nbins = 72;
1019 xmin = 0;
1020 xmax = 2*pi;
1021 break;
1022
1023 case 5:
1024 label = "Cluster Energy";
1025 nbins = 1000;
1026 xmin = 0;
1027 xmax = 10;
1028 break;
1029
1030 case 6:
1031 label = "Cluster Eta";
1032 nbins = 24;
1033 xmin = -1.2;
1034 xmax = 1.2;
1035 break;
1036
1037 case 7:
1038 label = "Cluster Phi";
1039 nbins = 72;
1040 xmin = 0;
1041 xmax = 2*pi;
1042 break;
1043
1044 case 8:
1045 label = "Is EMCalCluster";
1046 nbins = 2;
1047 xmin = 0;
1048 xmax = 2;
1049 break;
1050
1051 case 9:
1052 label = "Number of Matched Tracks to Cluster";
1053 nbins = 60;
1054 xmin = 0;
1055 xmax = 60;
1056 break;
1057
1058 case 10:
1059 label = "Track Pt";
1060 nbins = 750;
1061 xmin = 0;
1062 xmax = 75;
1063 break;
1064
1065 case 11:
1066 label = "Track Eta";
1067 nbins = 24;
1068 xmin = -1.2;
1069 xmax = 1.2;
1070 break;
1071
1072 case 12:
1073 label= "Track Phi";
1074 nbins = 72;
1075 xmin = 0;
1076 xmax = 2*pi;
1077 break;
1078
1079 case 13:
1080 label="Is Track EMCal";
1081 nbins = 2;
1082 xmin = 0;
1083 xmax = 2;
1084 break;
1085
1086 case 14:
1087 label = "Get Track EMCal Cluster";
1088 nbins = 100;
1089 xmin = 0;
1090 xmax = 100;
1091 break;
1092
1093 case 15:
1094 label = "Track Matched Phi";
1095 nbins = 72;
1096 xmin = 0;
1097 xmax = 2*pi;
1098
1099
1100
1101
1102 } // end of switch
1103} // end of getting dim-params
1104
5e5a2b89 1105
1106