[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderPlots.cxx


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


/*
 
$Log$
Revision 1.1.1.1  2003/05/29 18:56:58  horner
Initial import - Mark


*/

//_________________________________________________________________________
//  Class for Filling JetFinder Plots
//
//*-- Author: Mark Horner (LBL/UCT)
//
//


#include "TMath.h"
#include "AliEMCALJetFinderPlots.h"

ClassImp(AliEMCALJetFinderPlots)
	
AliEMCALJetFinderPlots::AliEMCALJetFinderPlots()
{
  fInitialised = kFALSE;
  fNominalEnergy = 0.0;
  fConeRadius = 0.5;
}

void AliEMCALJetFinderPlots::InitPlots()
{
	
    hFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
    hFragmFcn->Sumw2();
    hPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",100,0,1);
    hPartonFragmFcn->Sumw2();
    hPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
    hPartonJT->Sumw2();
    hPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
    hPartonPL->Sumw2();
    hJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
    hJetJT->Sumw2();
    hJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
    hJetPL->Sumw2();
    hJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
    hJetEt->Sumw2();
    hJetEta = new	TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
    hJetEta->Sumw2();
    hJetPhi = new 	TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
    hJetPhi->Sumw2();
    hPartonEta = new	TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
    hPartonEta->Sumw2();
    hPartonPhi = new 	TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
    hPartonPhi->Sumw2();
    hEtaDiff = new	TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
    hEtaDiff->Sumw2();
    hPhiDiff  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
    hPhiDiff->Sumw2();
    hNJets = new	TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
    hNJets->Sumw2();
    hEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
    hEtaPhiSpread->Sumw2();    
  hNJets->SetXTitle("N_{jets}^{reco}/event");
  hNJets->SetYTitle("N_{events}");

  //Jet properties
  hJetEt->SetFillColor(16);
  hJetEt->SetXTitle("E_{T}^{reco}");
  
  hJetEta->SetFillColor(16);
  hJetEta->SetXTitle("#eta_{jet}^{reco}");
  
  hJetPhi->SetFillColor(16);
  hJetPhi->SetXTitle("#phi_{jet}^{reco}");
  
  hPartonEta->SetFillColor(16);
  hPartonEta->SetXTitle("#eta_{parton}");
  
  hPartonPhi->SetFillColor(16);
  hPartonPhi->SetXTitle("#phi_{parton}");

  hPartonPL->SetXTitle("p (GeV/c)");
  hPartonJT->SetXTitle("p (GeV/c)");

  hPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");

  //Jet component properties

  hJetPL->SetXTitle("p (GeV/c)");
  hJetJT->SetXTitle("p (GeV/c)");
  hFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
  hPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");

  hEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
  hPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
  hEtaPhiSpread->SetXTitle("#eta");
  hEtaPhiSpread->SetYTitle("#phi");
  fInitialised = kTRUE;	
  
}

AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots(){

delete    hFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
delete   hPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
delete    hPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
delete    hJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
delete    hJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
delete    hJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
delete    hJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
delete    hJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
delete    hPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
delete    hPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
delete    hEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
delete    hPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
delete    hNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
delete 	  hEtaPhiSpread;	    
}	

void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output)
{
if (!fInitialised) InitPlots();	
  fOutput = output;
	
  Int_t nPartons = fOutput->GetNPartons();
  hNJets->Fill(fOutput->GetNJets());
  if (fOutput->GetNJets()!=1) return;
  AliEMCALParton* parton;
  AliEMCALJet* jet; 
  jet = fOutput->GetJet(0);
  
  hJetEt->Fill(jet->Energy());
  hJetEta->Fill(jet->Eta()  );
  hJetPhi->Fill(jet->Phi() );
  if (nPartons ==0) return;
  parton = fOutput->GetParton(0);
  hPartonEta->Fill( parton->Eta() );
  hPartonPhi->Fill( parton->Phi() );

  //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
  hEtaDiff->Fill( jet->Eta() - parton->Eta() );
  hPhiDiff->Fill( jet->Phi() - parton->Phi() );
  hEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
 /* 
  Float_t *pt,*phi,*eta;
  Int_t *pdg;
  pt  = new Float_t[parton->GetNTracks()];
  eta = new Float_t[parton->GetNTracks()];
  phi = new Float_t[parton->GetNTracks()];
  pdg = new Int_t[parton->GetNTracks()];*/


 Float_t pt[2000];
 Float_t eta[2000];
 Float_t phi[2000];
 Int_t pdg[2000];
  
  parton->GetTrackList(pt,eta,phi,pdg);
  for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
  {
      if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
           (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
      Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
      Double_t rt = 2.0*atan(exp(-parton->Eta()));
      Double_t ctt = cos(tt);
      Double_t crt = cos(rt);
      Double_t stt = sin(tt);
      Double_t srt = sin(rt);
      Double_t ctp = cos(phi[iT]);
      Double_t crp = cos(parton->Phi());
      Double_t stp = sin(phi[iT]);
      Double_t srp = sin(parton->Phi());
      Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
      Double_t correctp = pt[iT]/stt; 	
      hPartonPL->Fill( correctp*cos(alpha));
      if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
      (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
    	      hPartonJT->Fill( correctp*sin(alpha));
      if (fNominalEnergy == 0.0) {
	      hPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
      }else 
      {
              hPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
      }
  }// loop over tracks

/*
  pt  = new Float_t[jet->NTracks()];
  eta = new Float_t[jet->NTracks()];
  phi = new Float_t[jet->NTracks()];
  pdg = new Int_t[jet->NTracks()];*/
  jet->TrackList(pt,eta,phi,pdg);
  for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
  {
      	  Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
	  Double_t rt = 2.0*atan(exp(-jet->Eta()));
	  Double_t ctt = cos(tt);                   
	  Double_t crt = cos(rt);                         
	  Double_t stt = sin(tt);                               
	  Double_t srt = sin(rt);
	  Double_t ctp = cos(phi[iT]);
	  Double_t crp = cos(jet->Phi());
	  Double_t stp = sin(phi[iT]);
	  Double_t srp = sin(jet->Phi());
	  Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
	  Double_t correctp = pt[iT]/stt;
	  hJetPL->Fill( correctp*cos(alpha));      
	  if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +  
			  (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )   
		  hJetJT->Fill( correctp*sin(alpha));
	  if (fNominalEnergy==0.0){
		  hFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
	  } else
	  {
                  hFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
	  }
  }// loop over tracks
  
}
Commit	Line	Data
f7d5860b	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17
	18	/*
	19
	20	$Log$
	21	Revision 1.1.1.1 2003/05/29 18:56:58 horner
	22	Initial import - Mark
	23
	24
	25	*/
	26
	27	//_________________________________________________________________________
	28	// Class for Filling JetFinder Plots
	29	//
	30	//*-- Author: Mark Horner (LBL/UCT)
	31	//
	32	//
	33
	34
	35	#include "TMath.h"
	36	#include "AliEMCALJetFinderPlots.h"
	37
	38	ClassImp(AliEMCALJetFinderPlots)
	39
	40	AliEMCALJetFinderPlots::AliEMCALJetFinderPlots()
	41	{
	42	fInitialised = kFALSE;
	43	fNominalEnergy = 0.0;
	44	fConeRadius = 0.5;
	45	}
	46
	47	void AliEMCALJetFinderPlots::InitPlots()
	48	{
	49
	50	hFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
	51	hFragmFcn->Sumw2();
	52	hPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",100,0,1);
	53	hPartonFragmFcn->Sumw2();
	54	hPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
	55	hPartonJT->Sumw2();
	56	hPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
	57	hPartonPL->Sumw2();
	58	hJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
	59	hJetJT->Sumw2();
	60	hJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
	61	hJetPL->Sumw2();
	62	hJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
	63	hJetEt->Sumw2();
	64	hJetEta = new TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
65	hJetEta->Sumw2();
66	hJetPhi = new TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
67	hJetPhi->Sumw2();
68	hPartonEta = new TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
69	hPartonEta->Sumw2();
70	hPartonPhi = new TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
71	hPartonPhi->Sumw2();
72	hEtaDiff = new TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
73	hEtaDiff->Sumw2();
74	hPhiDiff = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
75	hPhiDiff->Sumw2();
76	hNJets = new TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
77	hNJets->Sumw2();
78	hEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
79	hEtaPhiSpread->Sumw2();
80	hNJets->SetXTitle("N_{jets}^{reco}/event");
81	hNJets->SetYTitle("N_{events}");
82
83	//Jet properties
84	hJetEt->SetFillColor(16);
85	hJetEt->SetXTitle("E_{T}^{reco}");
86
87	hJetEta->SetFillColor(16);
88	hJetEta->SetXTitle("#eta_{jet}^{reco}");
89
90	hJetPhi->SetFillColor(16);
91	hJetPhi->SetXTitle("#phi_{jet}^{reco}");
92
93	hPartonEta->SetFillColor(16);
94	hPartonEta->SetXTitle("#eta_{parton}");
95
96	hPartonPhi->SetFillColor(16);
97	hPartonPhi->SetXTitle("#phi_{parton}");
98
99	hPartonPL->SetXTitle("p (GeV/c)");
100	hPartonJT->SetXTitle("p (GeV/c)");
101
102	hPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");
103
104	//Jet component properties
105
106	hJetPL->SetXTitle("p (GeV/c)");
107	hJetJT->SetXTitle("p (GeV/c)");
108	hFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
109	hPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
110
111	hEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
112	hPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
113	hEtaPhiSpread->SetXTitle("#eta");
114	hEtaPhiSpread->SetYTitle("#phi");
115	fInitialised = kTRUE;
116
117	}
118
119	AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots(){
120
121	delete hFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
122	delete hPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
123	delete hPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
124	delete hJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
125	delete hJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
126	delete hJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
127	delete hJetEta;// = new TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
128	delete hJetPhi;// = new TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
129	delete hPartonEta;// = new TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
130	delete hPartonPhi;// = new TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
131	delete hEtaDiff;// = new TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
132	delete hPhiDiff;// = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
133	delete hNJets;// = new TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
134	delete hEtaPhiSpread;
135	}
136
137	void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output)
138	{
139	if (!fInitialised) InitPlots();
140	fOutput = output;
141
142	Int_t nPartons = fOutput->GetNPartons();
143	hNJets->Fill(fOutput->GetNJets());
144	if (fOutput->GetNJets()!=1) return;
145	AliEMCALParton* parton;
146	AliEMCALJet* jet;
147	jet = fOutput->GetJet(0);
148
149	hJetEt->Fill(jet->Energy());
150	hJetEta->Fill(jet->Eta() );
151	hJetPhi->Fill(jet->Phi() );
152	if (nPartons ==0) return;
153	parton = fOutput->GetParton(0);
154	hPartonEta->Fill( parton->Eta() );
155	hPartonPhi->Fill( parton->Phi() );
156
157	//hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
158	hEtaDiff->Fill( jet->Eta() - parton->Eta() );
159	hPhiDiff->Fill( jet->Phi() - parton->Phi() );
160	hEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
161	/*
162	Float_t pt,phi,*eta;
163	Int_t *pdg;
164	pt = new Float_t[parton->GetNTracks()];
165	eta = new Float_t[parton->GetNTracks()];
166	phi = new Float_t[parton->GetNTracks()];
167	pdg = new Int_t[parton->GetNTracks()];*/
168
169
170
171	Float_t pt[2000];
172	Float_t eta[2000];
173	Float_t phi[2000];
174	Int_t pdg[2000];
175
176	parton->GetTrackList(pt,eta,phi,pdg);
177	for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
178	{
179	if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
180	(phi[iT]-parton->Phi())(phi[iT]-parton->Phi()) >fConeRadius fConeRadius ) continue;
181	Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
182	Double_t rt = 2.0*atan(exp(-parton->Eta()));
183	Double_t ctt = cos(tt);
184	Double_t crt = cos(rt);
185	Double_t stt = sin(tt);
186	Double_t srt = sin(rt);
187	Double_t ctp = cos(phi[iT]);
188	Double_t crp = cos(parton->Phi());
189	Double_t stp = sin(phi[iT]);
190	Double_t srp = sin(parton->Phi());
191	Double_t alpha = acos(crpctpsrtstt+srpstpsrtstt+crt*ctt);
192	Double_t correctp = pt[iT]/stt;
193	hPartonPL->Fill( correctp*cos(alpha));
194	if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
195	(parton->Phi()-phi[iT])(parton->Phi()-phi[iT]) < 0.20.2 )
196	hPartonJT->Fill( correctp*sin(alpha));
197	if (fNominalEnergy == 0.0) {
198	hPartonFragmFcn->Fill( correctp*sin(tt)/parton->Energy() );
199	}else
200	{
201	hPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
202	}
203	}// loop over tracks
204
205	/*
206	pt = new Float_t[jet->NTracks()];
207	eta = new Float_t[jet->NTracks()];
208	phi = new Float_t[jet->NTracks()];
209	pdg = new Int_t[jet->NTracks()];*/
210	jet->TrackList(pt,eta,phi,pdg);
211	for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
212	{
213	Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
214	Double_t rt = 2.0*atan(exp(-jet->Eta()));
215	Double_t ctt = cos(tt);
216	Double_t crt = cos(rt);
217	Double_t stt = sin(tt);
218	Double_t srt = sin(rt);
219	Double_t ctp = cos(phi[iT]);
220	Double_t crp = cos(jet->Phi());
221	Double_t stp = sin(phi[iT]);
222	Double_t srp = sin(jet->Phi());
223	Double_t alpha = acos(crpctpsrtstt+srpstpsrtstt+crt*ctt);
224	Double_t correctp = pt[iT]/stt;
225	hJetPL->Fill( correctp*cos(alpha));
226	if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +
227	(jet->Phi()-phi[iT])(jet->Phi()-phi[iT]) < 0.20.2 )
228	hJetJT->Fill( correctp*sin(alpha));
229	if (fNominalEnergy==0.0){
230	hFragmFcn->Fill( correctp*sin(tt)/parton->Energy() );
231	} else
232	{
233	hFragmFcn->Fill( correctp*sin(tt)/fNominalEnergy );
234	}
235	}// loop over tracks
236
237	}
238
239
240