1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 // Jet fragmentation transverse momentum (j_T) analysis task
21 #include <TClonesArray.h>
27 #include <TLorentzVector.h>
29 #include <TGraphErrors.h>
34 #include "AliCentrality.h"
38 #include "AliVCluster.h"
39 #include "AliAODCaloCluster.h"
40 #include "AliESDCaloCluster.h"
41 #include "AliVTrack.h"
42 #include "AliEmcalJet.h"
43 #include "AliRhoParameter.h"
45 #include "AliJetContainer.h"
46 #include "AliParticleContainer.h"
47 #include "AliClusterContainer.h"
48 #include "AliPicoTrack.h"
50 #include "AliAnalysisTaskJetJTJT.h"
53 ClassImp(AliAnalysisTaskJetJTJT)
55 //________________________________________________________________________
56 AliAnalysisTaskJetJTJT::AliAnalysisTaskJetJTJT() :
57 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetJTJT", kTRUE),
63 fHistBackgroundDone(0),
74 fhTrackingEfficiency(0),
83 fTrackArrayName("nonejk"),
90 // Default constructor.
92 fHistTracksPt = new TH1*[fNcentBins];
93 fHistTracksJt = new TH1*[fNcentBins];
94 fHistClustersPt = new TH1*[fNcentBins];
95 fHistLeadingJetPt = new TH1*[fNcentBins];
96 fHistJetsPt = new TH1**[fNcentBins];
97 fHistBackgroundDone = new TH1**[fNcentBins];
98 fHistJTPta = new TH1***[fNcentBins];
99 fHistLogJTPta = new TH1***[fNcentBins];
100 fHistJTPta_all = new TH1***[fNcentBins];
101 fHistJTBg = new TH1***[fNcentBins];
102 fHistLogJTBg = new TH1***[fNcentBins];
103 fHistBgMulti = new TH1**[fNcentBins];
104 fHistBgPt = new TH1**[fNcentBins];
105 fHistJetEta = new TH1**[fNcentBins];
106 fHistJetMulti = new TH1**[fNcentBins];
107 fHistJetTracksPt = new TH1**[fNcentBins];
108 fhTrackingEfficiency = new TProfile*[fNcentBins];
109 //CentBinBorders = new Double_t[10];
112 for (Int_t i = 0; i < fNcentBins; i++) {
114 fHistLogJTPta[i] = 0;
115 fHistJTPta_all[i] = 0;
118 fHistBackgroundDone[i] = 0;
119 fHistTracksPt[i] = 0;
120 fHistTracksJt[i] = 0;
121 fHistClustersPt[i] = 0;
122 fHistLeadingJetPt[i] = 0;
127 fHistJetMulti[i] = 0;
128 fHistJetTracksPt[i] = 0;
129 fhTrackingEfficiency[i] = 0;
132 /*for(Int_t i = 0; i < nJetsConts; i++){
135 SetMakeGeneralHistograms(kTRUE);
138 //________________________________________________________________________
139 AliAnalysisTaskJetJTJT::AliAnalysisTaskJetJTJT(const char *name) :
140 AliAnalysisTaskEmcalJet(name, kTRUE),
144 fHistLeadingJetPt(0),
146 fHistBackgroundDone(0),
157 fhTrackingEfficiency(0),
164 fCaloClustersCont(0),
166 fTrackArrayName("nonejk"),
171 // Standard constructor.
172 fHistTracksPt = new TH1*[fNcentBins];
173 fHistTracksJt = new TH1*[fNcentBins];
174 fHistClustersPt = new TH1*[fNcentBins];
175 fHistLeadingJetPt = new TH1*[fNcentBins];
176 fHistJetsPt = new TH1**[fNcentBins];
177 fHistBackgroundDone = new TH1**[fNcentBins];
178 fHistJTPta = new TH1***[fNcentBins];
179 fHistLogJTPta = new TH1***[fNcentBins];
180 fHistJTPta_all = new TH1***[fNcentBins];
181 fHistJTBg = new TH1***[fNcentBins];
182 fHistLogJTBg = new TH1***[fNcentBins];
183 fHistBgMulti = new TH1**[fNcentBins];
184 fHistBgPt = new TH1**[fNcentBins];
185 fHistJetEta = new TH1**[fNcentBins];
186 fHistJetMulti = new TH1**[fNcentBins];
187 fHistJetTracksPt = new TH1**[fNcentBins];
188 fhTrackingEfficiency = new TProfile*[fNcentBins];
189 //CentBinBorders = new Double_t[10];
192 for (Int_t i = 0; i < fNcentBins; i++) {
194 fHistLogJTPta[i] = 0;
195 fHistJTPta_all[i] = 0;
198 fHistBackgroundDone[i] = 0;
199 fHistTracksPt[i] = 0;
200 fHistTracksJt[i] = 0;
201 fHistClustersPt[i] = 0;
202 fHistLeadingJetPt[i] = 0;
207 fHistJetMulti[i] = 0;
208 fHistJetTracksPt[i] = 0;
209 fhTrackingEfficiency[i] = 0;
212 /*for(Int_t i = 0; i < nJetsConts; i++){
215 SetMakeGeneralHistograms(kTRUE);
218 //________________________________________________________________________
219 AliAnalysisTaskJetJTJT::~AliAnalysisTaskJetJTJT()
225 void AliAnalysisTaskJetJTJT::setCentBinBorders( int n, Double_t *c){
228 cout << "AliAnalysisTaskJetJTJT::setCentBinBorders: " << endl;
230 for(int i= 0 ; i < fNcentBins; i++){
231 CentBinBorders[i]= c[i];
233 cout << CentBinBorders[i] << endl;
237 void AliAnalysisTaskJetJTJT::setTriggPtBorders( int n, Double_t *c){
240 cout << "AliAnalysisTaskJetJTJT::setTriggPtBorders: " << endl;
241 for(int i= 0 ; i < fNpttBins; i++){
242 TriggPtBorders[i]= c[i];
244 cout << TriggPtBorders[i] << endl;
248 void AliAnalysisTaskJetJTJT::setAssocPtBorders( int n, Double_t *c){
251 cout << "AliAnalysisTaskJetJTJT::setAssocPtBorders: " << endl;
252 for(int i= 0 ; i < fNptaBins; i++){
253 AssocPtBorders[i]= c[i];
255 cout << AssocPtBorders[i] << endl;
260 //________________________________________________________________________
261 void AliAnalysisTaskJetJTJT::UserCreateOutputObjects()
263 // Create user output.
265 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
267 cout << "Creating Histograms" << endl;
269 fJetsCont = GetJetContainer(0);
270 /*for(int i = 0; i < nJetsConts ; i++){
271 fJetsConts[0] = GetJetContainer(i);
273 /*if(fJetsCont) { //get particles and clusters connected to jets
274 fTracksCont = fJetsCont->GetParticleContainer();
275 fCaloClustersCont = fJetsCont->GetClusterContainer();
276 } else { //no jets, just analysis tracks and clusters
277 fTracksCont = GetParticleContainer(0);
278 fCaloClustersCont = GetClusterContainer(0);
280 fTracksCont = GetParticleContainer(0);
281 fCaloClustersCont = GetClusterContainer(0);
282 fTracksCont->SetClassName("AliVTrack");
283 fCaloClustersCont->SetClassName("AliAODCaloCluster");
286 //Int_t fMinBinJt = 0;
287 //Int_t fMaxBinJt = 5;
290 double LogBinsJt[NBINSJt+1], LimLJt=0.01, LimHJt=10;
291 double logBWJt = (TMath::Log(LimHJt)-TMath::Log(LimLJt))/(NBINSJt-1);
293 for(int ij=1;ij<=NBINSJt;ij++) LogBinsJt[ij]=LimLJt*exp(ij*logBWJt);
296 double LimLJtW=TMath::Log(0.01), LimHJtW=TMath::Log(10);
298 //==== Efficiency ====
300 cout << "AliAnalysisTaskJetJTJT::UserCreateOutputObjects: Creating efficiency" << endl;
301 fEfficiency = new JTJTEfficiency;
302 fEfficiency->SetMode( 1 ); // 0:NoEff, 1:Period 2:RunNum 3:Auto
303 fEfficiency->SetDataPath("alien:///alice/cern.ch/user/d/djkim/legotrain/efficieny/data"); // Efficiency root file location local or alien
307 for (Int_t ic = 0; ic < fNcentBins; ic++) {
308 if (fParticleCollArray.GetEntriesFast()>0) {
309 histname = "fHistTracksPt_";
311 fHistTracksPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt / 4);
312 fHistTracksPt[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
313 fHistTracksPt[ic]->GetYaxis()->SetTitle("counts");
314 fOutput->Add(fHistTracksPt[ic]);
317 if (fParticleCollArray.GetEntriesFast()>0) {
318 histname = "fHistTracksJt_";
320 fHistTracksJt[ic] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
321 fHistTracksJt[ic]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
322 fHistTracksJt[ic]->GetYaxis()->SetTitle("counts");
323 fOutput->Add(fHistTracksJt[ic]);
326 if (fParticleCollArray.GetEntriesFast()>0) {
327 histname = "fhTrackingEfficiency_";
329 fhTrackingEfficiency[ic] = new TProfile(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
330 fhTrackingEfficiency[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
331 fhTrackingEfficiency[ic]->GetYaxis()->SetTitle("counts");
332 fOutput->Add(fhTrackingEfficiency[ic]);
334 fHistJTPta[ic] = new TH1**[fNpttBins];
335 fHistLogJTPta[ic] = new TH1**[fNpttBins];
336 fHistJTPta_all[ic] = new TH1**[fNpttBins];
337 fHistJTBg[ic] = new TH1**[fNpttBins];
338 fHistLogJTBg[ic] = new TH1**[fNpttBins];
339 fHistJetsPt[ic] = new TH1*[fNpttBins];
340 fHistBackgroundDone[ic] = new TH1*[fNpttBins];
341 fHistBgMulti[ic] = new TH1*[fNpttBins];
342 fHistBgPt[ic] = new TH1*[fNpttBins];
343 fHistJetEta[ic] = new TH1*[fNpttBins];
344 fHistJetMulti[ic] = new TH1*[fNpttBins];
345 fHistJetTracksPt[ic] = new TH1*[fNpttBins];
346 for(Int_t j=0; j < fNpttBins; j++){
347 fHistJTPta[ic][j] = new TH1*[fNptaBins];
348 fHistLogJTPta[ic][j] = new TH1*[fNptaBins];
349 fHistJTPta_all[ic][j] = new TH1*[fNptaBins];
350 fHistJTBg[ic][j] = new TH1*[fNptaBins];
351 fHistLogJTBg[ic][j] = new TH1*[fNptaBins];
352 for(Int_t k=0; k < fNptaBins; k++){
353 fHistJTPta[ic][j][k] = 0;
354 fHistLogJTPta[ic][j][k] = 0;
355 fHistJTPta_all[ic][j][k] = 0;
356 fHistJTBg[ic][j][k] = 0;
357 fHistLogJTBg[ic][j][k] = 0;
359 fHistJetsPt[ic][j] = 0;
360 fHistBackgroundDone[ic][j] = 0;
361 fHistBgMulti[ic][j] = 0;
362 fHistBgPt[ic][j] = 0;
363 fHistJetEta[ic][j] = 0;
364 fHistJetMulti[ic][j] =0;
365 fHistJetTracksPt[ic][j] = 0;
369 if (fParticleCollArray.GetEntriesFast()>0) {
370 for(Int_t iptt = 0 ; iptt < fNpttBins; iptt++){
371 for(Int_t ipta = 0 ; ipta < fNptaBins; ipta++){
372 histname = "hJTPtaD00C";
374 histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
376 cout << histname << endl;
377 fHistJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
378 fHistJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
379 fHistJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
380 fOutput->Add(fHistJTPta[ic][iptt][ipta]);
382 histname = "hLogJTPtaD00C";
384 histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
386 cout << histname << endl;
387 fHistLogJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
388 fHistLogJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
389 fHistLogJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
390 fOutput->Add(fHistLogJTPta[ic][iptt][ipta]);
392 histname = "hJTPta_allD00C";
394 histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
396 cout << histname << endl;
397 fHistJTPta_all[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
398 fHistJTPta_all[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
399 fHistJTPta_all[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
400 fOutput->Add(fHistJTPta_all[ic][iptt][ipta]);
404 histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
406 cout << histname << endl;
407 fHistJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
408 fHistJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
409 fHistJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
410 fOutput->Add(fHistJTBg[ic][iptt][ipta]);
412 histname = "hLogJTBg";
414 histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
416 cout << histname << endl;
417 fHistLogJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
418 fHistLogJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
419 fHistLogJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
420 fOutput->Add(fHistLogJTBg[ic][iptt][ipta]);
427 if (fClusterCollArray.GetEntriesFast()>0) {
428 histname = "fHistClustersPt_";
430 fHistClustersPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
431 fHistClustersPt[ic]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
432 fHistClustersPt[ic]->GetYaxis()->SetTitle("counts");
433 fOutput->Add(fHistClustersPt[ic]);
436 if (fJetCollArray.GetEntriesFast()>0) {
437 histname = "fHistLeadingJetPt_";
439 fHistLeadingJetPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
440 fHistLeadingJetPt[ic]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
441 fHistLeadingJetPt[ic]->GetYaxis()->SetTitle("counts");
442 fOutput->Add(fHistLeadingJetPt[ic]);
444 for(Int_t iptt = 0 ; iptt < fNpttBins; iptt++){
446 histname = "fHistJetsPt_";
447 histname += Form("C%02dT%02d", ic, iptt);
448 fHistJetsPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
449 fHistJetsPt[ic][iptt]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
450 fHistJetsPt[ic][iptt]->GetYaxis()->SetTitle("counts");
451 fOutput->Add(fHistJetsPt[ic][iptt]);
453 histname = "fHistBackgroundDone_";
454 histname += Form("C%02dT%02d", ic, iptt);;
455 fHistBackgroundDone[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 2, -1, 2);
456 fHistBackgroundDone[ic][iptt]->GetXaxis()->SetTitle("Number of jets");
457 fHistBackgroundDone[ic][iptt]->GetYaxis()->SetTitle("0 = not done, 1 = done");
458 fOutput->Add(fHistBackgroundDone[ic][iptt]);
460 histname = "fHistJetEta_";
461 histname += Form("C%02dT%02d", ic, iptt);
462 fHistJetEta[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, -5, 5);
463 fHistJetEta[ic][iptt]->GetXaxis()->SetTitle("#eta");
464 fHistJetEta[ic][iptt]->GetYaxis()->SetTitle("jets");
465 fOutput->Add(fHistJetEta[ic][iptt]);
467 histname = "fHistJetMulti_";
468 histname += Form("C%02dT%02d", ic, iptt);
469 fHistJetMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
470 fHistJetMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
471 fHistJetMulti[ic][iptt]->GetYaxis()->SetTitle("jets");
472 fOutput->Add(fHistJetMulti[ic][iptt]);
474 histname = "fHistBgMulti_";
475 histname += Form("C%02dT%02d", ic, iptt);
476 fHistBgMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
477 fHistBgMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
478 fHistBgMulti[ic][iptt]->GetYaxis()->SetTitle("Events");
479 fOutput->Add(fHistBgMulti[ic][iptt]);
481 histname = "fHistBgPt_";
482 histname += Form("C%02dT%02d", ic, iptt);
483 fHistBgPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins,fMinBinPt, fMaxBinPt/4);
484 fHistBgPt[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
485 fHistBgPt[ic][iptt]->GetYaxis()->SetTitle("jets");
486 fOutput->Add(fHistBgPt[ic][iptt]);
488 histname = "fHistJetTracksPt_";
490 histname += Form("C%02dT%02d", ic, iptt);
492 cout << histname << endl;
493 fHistJetTracksPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt/4);
494 fHistJetTracksPt[ic][iptt]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
495 fHistJetTracksPt[ic][iptt]->GetYaxis()->SetTitle("counts");
496 fOutput->Add(fHistJetTracksPt[ic][iptt]);
500 if (!(GetJetContainer()->GetRhoName().IsNull())) {
501 histname = "fHistJetsCorrPtArea_";
503 fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
504 fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
505 fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
506 fOutput->Add(fHistJetsCorrPtArea[i]);
512 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
515 //________________________________________________________________________
516 Bool_t AliAnalysisTaskJetJTJT::FillHistograms()
520 AliCentrality *aliCent = InputEvent()->GetCentrality();
523 //fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
524 fCent = aliCent->GetCentralityPercentile("V0M");
526 cout << "Centrality " << fCent << endl;
528 for(int ic = 0; ic < fNcentBins; ic++){
530 cout << "ic: " << ic << " / " << fNcentBins << endl;
531 cout << "Centrality bin " << fCentBin << endl;
532 cout << "Border: " << CentBinBorders[ic] << endl;
534 if(fCent > CentBinBorders[ic]){
538 //cout << "Centrality bin: " << fCentBin << endl;
540 AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
543 int fHadronSelectionCut = 5; //5=Hybrid cut
545 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
547 double ptt = track->Pt();
549 //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
550 //AliJBaseTrack *triggTr = (AliJBaseTrack*)fchargedHadronList->At(i);
552 double effCorr = 1./fEfficiency->GetCorrection(ptt, fHadronSelectionCut, fCent); // here you generate warning if ptt>30
553 //double effCorr = 1.;
554 fhTrackingEfficiency[fCentBin]->Fill( ptt, 1./effCorr );
555 //triggTr->SetTrackEff( 1./effCorr );
556 //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
558 if(ptt > 0 && 1.0/ptt > 0){
559 fHistTracksPt[fCentBin]->Fill(ptt,1./ptt*effCorr);
563 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
567 if (fCaloClustersCont) {
568 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
570 TLorentzVector nPart;
571 cluster->GetMomentum(nPart, fVertex);
572 fHistClustersPt[fCentBin]->Fill(nPart.Pt());
574 cluster = fCaloClustersCont->GetNextAcceptCluster();
578 Int_t fPttBin, fPtaBin;
583 //Int_t Njets = fJetsCont->GetNJets();
586 cout << "Number of Jets: " << Njets << endl;
589 //Make a array to hold jets to be tested in background jT
590 Float_t jetPhis[200] = {};
591 Float_t jetEtas[200] = {};
592 AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
595 //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
596 if(jet->Pt() > 5){ //Only consider jets with pT > 5 GeV
597 jetPhis[ij] = jet->Phi();
598 jetEtas[ij] = jet->Eta();
602 cout << "i: " << ij << " jetPhi: " << jetPhis[ij] << " jetEta: " << jetEtas[ij] << endl;
608 cout << "jetPt: " << jet->Pt() << " jetPhi: " << jet->Phi() << " jetEta: " << jet->Eta() << endl;
611 jet = fJetsCont->GetNextAcceptJet();
615 //fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
616 jet = fJetsCont->GetNextAcceptJet(0);
619 if(jet->Eta() < -0.4 || jet->Eta() > 0.4){
621 cout << "Jet outside eta range, Eta: " << jet->Eta() << endl;
622 jet = fJetsCont->GetNextAcceptJet();
625 //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
626 //Get the trigger pT bin
628 for(int iptt = 0 ; iptt < fNpttBins; iptt++){
629 if(jet->Pt() > TriggPtBorders[iptt]){
634 fHistJetEta[fCentBin][fPttBin]->Fill(jet->Eta());
635 if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
636 fHistJetsPt[fCentBin][fPttBin]->Fill(jet->Pt(),1.0/jet->Pt()); //Fill jet dN/(pT dpT)
639 /*if (fHistJetsCorrPtArea[fCentBin]) {
640 Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
641 fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
643 //Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
646 Int_t nTrack = jet->GetNumberOfTracks();
648 cout << "Number of tracks " << nTrack << " Jet Pt: " << jet->Pt() << endl;
649 fHistJetMulti[fCentBin][fPttBin]->Fill(nTrack);
650 for(Int_t it = 0; it < nTrack; it++ ){
651 AliVParticle *track = (AliVParticle*)jet->TrackAt( it, fTracks );
653 cout << "No Track found" << endl;
656 fPtaBin = 0; //Get the associated pT bin
657 for(int ipta = 0 ; ipta < fNptaBins; ipta++){
658 if(track->Pt() > AssocPtBorders[ipta]){
662 fHistJetTracksPt[fCentBin][fPttBin]->Fill(track->Pt());
664 cout << "Filling fHistJetTracksPt C" << fCentBin << " T" << fPttBin << endl;
665 /*Float_t dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz(); // p_track dot p_jet
666 Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz()); // track pT norm
667 Float_t normproduct = constp*jetp; // jet pT norm times track pT norm
668 Float_t costheta2 = dotproduct/normproduct;
669 //Float_t sintheta = sqrt(1-costheta2*costheta2);
670 Float_t jt = constp*sqrt(1-costheta2*costheta2);*/
671 Float_t jt = getJt(track,jet,0);
672 double effCorr = 1./fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent); // here you generate warning if ptt>30
673 //double effCorr = 1.;
674 if(jt > 0 && 1.0/jt > 0){
675 fHistTracksJt[fCentBin]->Fill(jt,1.0/jt*effCorr); //Fill dN/(djT jT)
676 fHistJTPta[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr); //Fill dN/(djT jT)
677 fHistLogJTPta[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0/jt*effCorr); //Fill logarithmic dN/(djT jT)
680 cout << "Filling JT C" << fCentBin << "T" << fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << 1.0/jt*effCorr << endl;
683 //Get Jet azimuth and rapidity of jet
684 Float_t jetAngle = jet->Phi();
685 Float_t jetRap = jet->Eta();
687 //Rotate jet angle for background cone
688 Float_t rotatedAngle = jetAngle+TMath::Pi()/2;
689 if(rotatedAngle > TMath::Pi()*2){
690 rotatedAngle = rotatedAngle- TMath::Pi()*2;
692 Float_t jetArea = jet->Area();
693 Float_t testRadius = TMath::Sqrt(jetArea/TMath::Pi());
697 //Test if there are jets in the background test cone
698 for(int i_j = 0; i_j < Njets; i_j++){
699 Float_t diffR = TMath::Sqrt(TMath::Power(jetPhis[i_j]-rotatedAngle,2)+TMath::Power(jetEtas[i_j]-jetRap,2));
701 cout << "i_j: " << i_j << " JetPhi: " << jetPhis[i_j] << " jetEta: " << jetEtas[i_j] << endl;
702 cout << "DiffR: " << diffR << " doBG: " << doBg <<endl;
704 if(diffR < testRadius *2){ //Jets muts be at least 2*cone radius away from the background cone axis
711 // Do jT for all particles in respect to jet axis
714 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
716 /*Float_t dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
717 Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
718 Float_t normproduct = constp*jetp;
719 Float_t costheta2 = dotproduct/normproduct;
720 //Float_t sintheta = sqrt(1-costheta2*costheta2);
721 Float_t jt = constp*sqrt(1-costheta2*costheta2);*/
722 Double_t jt = getJt(track,jet,0);
723 double effCorr = 1./fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent); // here you generate warning if ptt>30
724 if(jt > 0 && 1.0/jt > 0){
725 fHistJTPta_all[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr);
727 for(int ipta = 0 ; ipta < fNptaBins; ipta++){
728 if(track->Pt() > AssocPtBorders[ipta]){
733 //If background is to be filled
735 Float_t diffR = TMath::Sqrt(TMath::Power(track->Phi()-rotatedAngle,2)+TMath::Power(track->Eta()-jetRap,2));
736 //Particles in the rotated cone
737 if(diffR < testRadius){
739 fHistBgPt[fCentBin][fPttBin]->Fill(track->Pt());
740 /*dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
741 constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
742 normproduct = constp*jetp;
743 costheta2 = dotproduct/normproduct;
744 //sintheta = sqrt(1-costheta2*costheta2);
745 jt = constp*sqrt(1-costheta2*costheta2);*/
746 jt = getJt(track,jet,1);
747 if(jt > 0 && 1.0/jt > 0){
748 fHistJTBg[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr);
749 fHistLogJTBg[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0/jt*effCorr);
752 cout << "Filling Background C" << fCentBin << "T" << fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << 1.0/jt*effCorr << endl;
757 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
760 fHistBgMulti[fCentBin][fPttBin]->Fill(counter);
764 fHistBackgroundDone[fCentBin][fPttBin]->Fill(1);
766 fHistBackgroundDone[fCentBin][fPttBin]->Fill(0);
770 jet = fJetsCont->GetNextAcceptJet();
773 jet = fJetsCont->GetLeadingJet();
775 if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
776 fHistLeadingJetPt[fCentBin]->Fill(jet->Pt(),1.0/jet->Pt());
781 CheckClusTrackMatching();
787 //-----------------------------------------------------------------------
788 Double_t AliAnalysisTaskJetJTJT::getJt(AliVTrack *track, AliEmcalJet *jet,int reverse){
789 Float_t dotproduct = 0;
790 Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
792 dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
794 dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
796 Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
797 Float_t normproduct = constp*jetp;
798 Float_t costheta2 = dotproduct/normproduct;
799 //Float_t sintheta = sqrt(1-costheta2*costheta2);
800 Float_t jt = constp*sqrt(1-costheta2*costheta2);
804 Double_t AliAnalysisTaskJetJTJT::getJt(AliVParticle *track, AliEmcalJet *jet,int reverse){
805 Float_t dotproduct = 0;
806 Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
808 dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
810 dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
812 Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
813 Float_t normproduct = constp*jetp;
814 Float_t costheta2 = dotproduct/normproduct;
815 //Float_t sintheta = sqrt(1-costheta2*costheta2);
816 Float_t jt = constp*sqrt(1-costheta2*costheta2);
820 //________________________________________________________________________
821 void AliAnalysisTaskJetJTJT::CheckClusTrackMatching()
824 if(!fTracksCont || !fCaloClustersCont)
830 //Get closest cluster to track
831 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
833 //Get matched cluster
834 Int_t emc1 = track->GetEMCALcluster();
835 if(fCaloClustersCont && emc1>=0) {
836 AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
838 AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
839 //fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
842 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
845 //Get closest track to cluster
846 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
848 TLorentzVector nPart;
849 cluster->GetMomentum(nPart, fVertex);
850 fHistClustersPt[fCentBin]->Fill(nPart.Pt());
853 AliVTrack *mt = NULL;
854 AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
856 if(acl->GetNTracksMatched()>1)
857 mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
860 AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
861 Int_t im = ecl->GetTrackMatchedIndex();
862 if(fTracksCont && im>=0) {
863 mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
867 AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
868 //fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
873 Int_t emc1 = mt->GetEMCALcluster();
874 Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
875 AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
876 AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
877 Printf("deta: %f dphi: %f",deta,dphi);
881 cluster = fCaloClustersCont->GetNextAcceptCluster();
885 //________________________________________________________________________
886 /*AliJetContainer* AliAnalysisTaskJetJTJT::AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius) {
888 AliAnalysisTaskEmcalJet::ExecOnce();
890 AliJetContainer *cont = 0x0;
891 cont = AliAnalysisTaskEmcalJet::AddJetContainer(n,defaultCutType,jetRadius);
896 //________________________________________________________________________
897 void AliAnalysisTaskJetJTJT::ExecOnce() {
900 cout << "AliAnalysisTaskJetJTJT::ExecOnce(): " << endl;
901 cout << "Get fTracks from " << fTrackArrayName.Data() << endl;
903 fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
905 AliAnalysisTaskEmcalJet::ExecOnce();
907 cout << "Efficiency: Set Run Period Name " << runPeriod << endl;
908 fEfficiency->SetPeriodName(runPeriod);
910 cout << "Efficiency: Set Run number " << InputEvent()->GetRunNumber() << endl;
911 fEfficiency->SetRunNumber( InputEvent()->GetRunNumber() ); //TODO Get run Number
913 cout << "Efficiency: Load()" << endl;
916 cout << "fEfficiency loaded" << endl;
920 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
921 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
922 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
926 //________________________________________________________________________
927 Bool_t AliAnalysisTaskJetJTJT::Run()
929 // Run analysis code here, if needed. It will be executed before FillHistograms().
931 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
934 //________________________________________________________________________
935 void AliAnalysisTaskJetJTJT::Terminate(Option_t *)
937 // Called once at the end of the analysis.
942 //________________________________________________________________________
943 JTJTEfficiency::JTJTEfficiency():
956 for (Int_t i=0; i < 3; i++) fEffDir[i] = 0;
960 JTJTEfficiency::JTJTEfficiency(const JTJTEfficiency& obj) :
962 fPeriod(obj.fPeriod),
963 fDataPath(obj.fDataPath),
965 fPeriodStr(obj.fPeriodStr),
966 fMCPeriodStr(obj.fMCPeriodStr),
967 fRunNumber(obj.fRunNumber),
969 fInputRootName(obj.fInputRootName),
970 fInputRoot(obj.fInputRoot),
971 fCentBinAxis(obj.fCentBinAxis)
973 // copy constructor TODO: handling of pointer members
975 for (Int_t i=0; i < 3; i++) fEffDir[i] = obj.fEffDir[i];
978 JTJTEfficiency& JTJTEfficiency::operator=(const JTJTEfficiency& obj){
979 // equal sign operator TODO: content
985 //________________________________________________________________________
986 double JTJTEfficiency::GetCorrection( double pt, int icut , double cent ) const {
987 if( fMode == kNotUse ) return 1;
988 int icent = fCentBinAxis->FindBin( cent ) -1 ;
989 if( icent < 0 || icent > fCentBinAxis->GetNbins()-1 ) {
990 cout<<"J_WARNING : Centrality "<<cent<<" is out of CentBinBorder"<<endl;
993 // TODO error check for icent;
995 if( ! fCorrection[ivtx][icent][icut] ) {
996 cout<<"J_WARNING : No Eff Info "<<pt<<"\t"<<icut<<"\t"<<cent<<"\t"<<icent<<endl;
999 TGraphErrors * gr = fCorrection[ivtx][icent][icut];
1000 //=== TEMPERORY SETTING. IT will be removed soon.
1001 if( pt > 30 ) pt = 30; // Getting eff of 30GeV for lager pt
1002 double cor = gr->Eval(pt);
1003 if ( cor < 0.2 ) cor = 0.2;
1008 TString JTJTEfficiency::GetEffName() {
1010 1. kNotUse : no Load, efficiency is 1 always
1011 2. has fInputRootName : Load that or crash
1012 3. has fName : Load fName [+runnumber] or crash
1013 4. has runnumber : Find Good MC period from AliJRunTable, or crash
1014 3. has period : Find Good MC period from AliJRunTable, or crash
1018 if(fPeriodStr == "LHC10b"){
1019 fInputRootName = "Eff--LHC10b-LHC10d1-0-.root";
1021 if(fPeriodStr == "LHC10c"){
1022 fInputRootName = "Eff--LHC10c-LHC10d4-0-.root";
1024 if(fPeriodStr == "LHC10d"){
1025 fInputRootName = "Eff--LHC10d-LHC10f6a-0-.root";
1027 if(fPeriodStr == "LHC10e"){
1028 fInputRootName = "Eff--LHC10e-LHC10e20-0-.root";
1030 if(fPeriodStr == "LHC10h"){
1031 fInputRootName = "Eff--LHC10h-LHC11a10a_bis-0-.root";
1033 if(fPeriodStr == "LHC11a"){
1034 fInputRootName = "Eff--LHC11a-LHC11b10a-0-.root";
1036 if(fPeriodStr == "LHC13b"){
1037 fInputRootName = "Eff--LHC13b-LHC13b2-efix_p1-0-.root";
1040 if(fPeriodStr == "LHC13c"){
1041 fInputRootName = "Eff--LHC13c-LHC13b2-efix_p1-0-.root";
1043 if(fPeriodStr == "LHC13d"){
1044 fInputRootName = "Eff--LHC13d-LHC13b2-efix_p1-0-.root";
1046 if(fPeriodStr == "LHC13e"){
1047 fInputRootName = "Eff--LHC13e-LHC13b2-efix_p1-0-.root";
1050 return fInputRootName;
1053 TString JTJTEfficiency::GetEffFullName() {
1055 fInputRootName = fDataPath + "/" + fInputRootName;
1056 return fInputRootName;
1060 //________________________________________________________________________
1061 bool JTJTEfficiency::Load(){
1062 // Load Efficiency File based on fMode
1063 if( fMode == kNotUse ) {
1064 cout<<"J_WARNING : Eff Mode is \"NOTUSE\". eff is 1 !!!"<<endl;
1068 if (TString(fInputRootName).BeginsWith("alien:")) TGrid::Connect("alien:");
1069 fInputRoot = TFile::Open( fInputRootName);
1070 //fInputRoot = new TFile( fInputRootName,"READ");
1072 cout << "J_ERROR : %s does not exist" << fInputRootName << endl;
1076 //fEffDir[0] = (TDirectory*)fInputRoot->Get("EffRE");
1077 ///fEffDir[1] = (TDirectory*)fInputRoot->Get("EffMC");
1078 fEffDir[2] = (TDirectory*)fInputRoot->Get("Efficiency");
1079 //iif( fEffDir[0] && fEffDir[1] && fEffDir[2] )
1082 cout << "J_ERROR : Directory EFF is not exist"<<endl;
1086 fCentBinAxis = (TAxis*)fEffDir[2]->Get("CentralityBin");
1087 if( !fCentBinAxis ){
1088 cout << "J_ERROR : No CentralityBin in directory" << endl;
1094 int nCentBin = fCentBinAxis->GetNbins();
1095 for( int ivtx=0;ivtx<nVtx;ivtx++ ){
1096 for( int icent=0;icent<nCentBin;icent++ ){
1097 for( int icut=0;icut<kJNTrackCuts;icut++ ){
1098 fCorrection[ivtx][icent][icut]
1099 = (TGraphErrors*) fEffDir[2]->Get(Form("gCor%02d%02d%02d", ivtx,icent,icut));
1100 //cout<<"J_LOG : Eff graph - "<<Form("gCor%02d%02d%02d", ivtx,icent,icut)<<" - "<<g<<endl;
1104 cout<<"J_LOG : Eff file is "<<fInputRootName<<endl;
1105 cout<<"J_LOG : Eff Cent Bins are ";
1106 for( int i=0;i<=nCentBin;i++ ){
1107 //cout<<fCentBinAxis->GetXbins()->At(i)<<" ";