StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHiSpectra.cxx
1 /***************************************************************************
2  *
3  * $Id: StHiSpectra.cxx,v 1.4 2003/04/30 20:37:30 perev Exp $
4  *
5  * Author: Bum Choi, UT Austin, Apr 2002
6  *
7  ***************************************************************************
8  *
9  * Description: Class for making highpt inclusive spectra from highpt
10  * uDST's
11  *
12  ***************************************************************************
13  *
14  * $Log: StHiSpectra.cxx,v $
15  * Revision 1.4 2003/04/30 20:37:30 perev
16  * Warnings cleanup. Modified lines marked VP
17  *
18  * Revision 1.3 2002/05/31 21:58:29 jklay
19  * Updated analysis code to use new cut class
20  *
21  * Revision 1.2 2002/04/03 00:23:27 jklay
22  * Fixed private member access bugs in analysis code
23  *
24  * Revision 1.1 2002/04/02 20:05:18 jklay
25  * Bums analysis tools for highpt uDSTs
26  *
27  *
28  **************************************************************************/
29 #include "StHiSpectra.h"
30 
31 
32 const char* border = "*************************************************";
33 
34 //__________________
35 
36 StHiSpectra::StHiSpectra(const char* inputDir,
37  const char* outRootName,
38  const char* f)
39  : StHiBaseAnalysis(inputDir,outRootName)
40 
41 {
42 
43  mEfficiencyString = ""; mEfficiencyFileName = f;
44  mEfficiencyMap = 0;
45 }
46 
47 //__________________
48 
49 StHiSpectra::~StHiSpectra()
50 {
51 }
52 
53 //__________________
54 //
55 // reads in all the files in the directory
56 // by default, looks for files ending in minimc.root
57 //
58 
59 Int_t
60 StHiSpectra::initMore()
61 {
62 
63  // char name[100];
64  Int_t stat=0;
65 
66  TString fName = "fEfficiency";
67  TFile* file=0;
68  if(mEfficiencyFileName==""){
69  cout << "no efficiency file set " << endl;
70  }
71  else{
72  file=new TFile(mEfficiencyFileName.Data());
73  if(!file ||!file->IsOpen()){
74  cout << "Cannot find efficiency file : " << mEfficiencyFileName << endl;
75  }
76  else{
77  mEfficiencyMap = (TF2*)file->Get(fName.Data());
78  delete file;
79  }
80  }
81  if(!mEfficiencyMap){
82  cout << "Cannot find efficency function : " << fName.Data() <<endl;
83  cout << "WARNING-USING DUMMY EFFIENCY" << endl;
84  mEfficiencyMap = new TF2("fEfficiency","1");
85  }
86 
87 
88  cout << "\tefficiency file : " << mEfficiencyFileName << endl;
89  cout << "\tefficiency fcn title : " << mEfficiencyMap->GetTitle() << endl;
90 
91  /*
92  if(mEfficiencyString==""){
93  cout << "no efficiency string" << endl;
94  stat++;
95  }
96 
97  else{
98  // create the efficiency function
99  // x is eta, y is pt
100  cout << "efficiency : " << mEfficiencyString << endl;
101 
102  mEfficiencyMap = new TF2("fEfficiencyMap",mEfficiencyString,-.55,.55,1.6,6);
103  }
104  */
105  return stat;
106 }
107 
108 
109 void
110 StHiSpectra::initHistograms()
111 {
112  cout << "StHiSpectra::initHistograms()" << endl;
113 
114  //*********************
115 
116  using namespace Bin;
117 
118  //*********************
119 
120  // gStyle->SetPalette(1,0);
121 
122  char name[500];
123 
124  h1CentralityCut = new TH1D("h1CentralityCut","centrality cut",
125  nFlowCentBin,flowCentMin,flowCentMax);
126 
127  h1Centrality = new TH1D("h1Centrality","centrality",
128  nFlowCentBin,flowCentMin,flowCentMax);
129 
130  h1VertexZ = new TH1D("h1VertexZ","vertex z",
131  nVertexZEvtBin,vertexZEvtMin,vertexZEvtMax);
132 
133 // sprintf(title,"vertex z (%d<=zdc cent<=%d)",
134 // CutRc::mZdcCtbCent[1],CutRc::mZdcCtbCent[0]);
135 // h1VertexZCut = new TH1D("h1VertexZCut",title,
136 // nVertexZEvtBin,vertexZEvtMin,vertexZEvtMax);
137 
138  h1VertexZThin = new TH1D("h1VertexZThin","vertex z",
139  nVertexZEvtThinBin,vertexZEvtMin,vertexZEvtMax);
140 
141 // sprintf(title,"vertex r %d<=zdc cent<=%d)",
142 // CutRc::mZdcCtbCent[1],CutRc::mZdcCtbCent[0]);
143 // h1VertexRCut = new TH1D("h1VertexRCut",title,
144 // nVertexREvtBin,vertexREvtMin,vertexREvtMax);
145 
146  cout << "using bins" << endl;
147  for(int iBin=0; iBin<2; iBin++){
148  varBin[iBin].ptBinsAry = new TArrayD;
149  initPtAry(varBin[iBin].ptBinsAry,iBin);
150  int nBin = varBin[iBin].ptBinsAry->GetSize()-1;
151  varBin[iBin].nPtBinAry = nBin;;
152  cout << "Bin " << iBin << " : " << endl;
153  for(int i=0; i<nBin; i++){
154  cout << varBin[iBin].ptBinsAry->At(i) << ", ";
155  }
156  cout << varBin[iBin].ptBinsAry->At(nBin) << endl;
157  }
158  //------------------------------------------------------------
159 
160  TString sPM[2] = { "Plus","Minus"};
161 
162  //********************************************************
163  // number of events used to scale the spectra
164 
165  h1NEvent = new TH1D("h1NEvent","h1NEvent",1,1,2);
166 
167  // eta cut
168 
169  h1EtaCut = new TH1D("h1EtaCut","h1EtaCut",2,0,2);
170 
171 
172  //********************************************************
173  // yield for both signs
174 
175  h2Yield
176  = new TH2D("h2Yield","h2Yield",
177  nEtaSmallBin,etaSmallMin,etaSmallMax,
178  nPtBin,ptMin,ptMax);
179  h2Yield->Sumw2();
180 
181 
182 
183 
184  //********************************************************
185  // var bin
186 
187  for(Int_t iBin=0; iBin<mNVarBin; iBin++){
188 
189  //**** both charge signs
190 
191  setName(name,"h1OneOverPt",iBin);
192  varBin[iBin].mean.h1OneOverPt =
193  new TH1D(name,name,varBin[iBin].nPtBinAry,
194  varBin[iBin].ptBinsAry->GetArray());
195  varBin[iBin].mean.h1OneOverPt->Sumw2();
196 
197  setName(name,"h1WeightedMean",iBin);
198  varBin[iBin].mean.h1WeightedMean =
199  new TH1D(name,name,varBin[iBin].nPtBinAry,
200  varBin[iBin].ptBinsAry->GetArray());
201  varBin[iBin].mean.h1WeightedMean->Sumw2();
202 
203 
204  //
205  // raw
206  //
207  setName(name,"h1Raw",iBin);
208  varBin[iBin].spec.h1Raw =
209  new TH1D(name,name,varBin[iBin].nPtBinAry,
210  varBin[iBin].ptBinsAry->GetArray());
211  varBin[iBin].spec.h1Raw->Sumw2();
212 
213  //
214  // efficiency corrected
215  //
216  setName(name,"h1EffCorrected",iBin);
217  varBin[iBin].spec.h1EffCorrected =
218  new TH1D(name,name,varBin[iBin].nPtBinAry,
219  varBin[iBin].ptBinsAry->GetArray());
220  varBin[iBin].spec.h1EffCorrected->Sumw2();
221 
222  //
223  // corrected (background + momentum resolution)
224  //
225  setName(name,"h1Corrected",iBin);
226  varBin[iBin].spec.h1Corrected =
227  new TH1D(name,name,varBin[iBin].nPtBinAry,
228  varBin[iBin].ptBinsAry->GetArray());
229  varBin[iBin].spec.h1Corrected->Sumw2();
230 
231 
232  //**** spectra plus, minus
233 
234  for(Int_t iCharge=0; iCharge<2; iCharge++){
235 
236  // mean stuff
237  //
238 
239  setName(name,"h1OneOverPt",iBin,sPM[iCharge].Data());
240  varBin[iBin].meanPM[iCharge].h1OneOverPt =
241  new TH1D(name,name,varBin[iBin].nPtBinAry,
242  varBin[iBin].ptBinsAry->GetArray());
243  varBin[iBin].meanPM[iCharge].h1OneOverPt->Sumw2();
244 
245  setName(name,"h1WeightedMean",iBin,sPM[iCharge].Data());
246  varBin[iBin].meanPM[iCharge].h1WeightedMean =
247  new TH1D(name,name,varBin[iBin].nPtBinAry,
248  varBin[iBin].ptBinsAry->GetArray());
249  varBin[iBin].meanPM[iCharge].h1WeightedMean->Sumw2();
250 
251  //
252  // raw
253  //
254  setName(name,"h1Raw",iBin,sPM[iCharge].Data());
255  varBin[iBin].specPM[iCharge].h1Raw =
256  new TH1D(name,name,varBin[iBin].nPtBinAry,
257  varBin[iBin].ptBinsAry->GetArray());
258  varBin[iBin].specPM[iCharge].h1Raw->Sumw2();
259 
260  //
261  // efficiency corrected
262  //
263  setName(name,"h1EffCorrected",iBin,sPM[iCharge].Data());
264  varBin[iBin].specPM[iCharge].h1EffCorrected =
265  new TH1D(name,name,varBin[iBin].nPtBinAry,
266  varBin[iBin].ptBinsAry->GetArray());
267  varBin[iBin].specPM[iCharge].h1EffCorrected->Sumw2();
268 
269  //
270  // corrected (background + momentum resolution)
271  //
272 
273  setName(name,"h1Corrected",iBin,sPM[iCharge].Data());
274  varBin[iBin].specPM[iCharge].h1Corrected =
275  new TH1D(name,name,varBin[iBin].nPtBinAry,
276  varBin[iBin].ptBinsAry->GetArray());
277  varBin[iBin].specPM[iCharge].h1Corrected->Sumw2();
278 
279 
280  }
281 
282  //***** ratios
283 
284  setName(name,"h1RawRatio",iBin);
285  varBin[iBin].h1RawRatio
286  = new TH1D(name,name,varBin[iBin].nPtBinAry,
287  varBin[iBin].ptBinsAry->GetArray());
288  varBin[iBin].h1RawRatio->Sumw2();
289 
290  setName(name,"h1CorrectedRatio",iBin);
291  varBin[iBin].h1CorrectedRatio
292  = new TH1D(name,name,varBin[iBin].nPtBinAry,
293  varBin[iBin].ptBinsAry->GetArray());
294  varBin[iBin].h1CorrectedRatio->Sumw2();
295 
296  } // varBin
297 }
298 
299 //_____________________
300 
301 void
302 StHiSpectra::trackLoop()
303 {
304  if(mDebug)
305  cout << "StHiSpectra::spectraTrackLoop()" << endl;
306 
307  Int_t nTrack = mHiMicroEvent->NTrack();
308  StHiMicroTrack* track;
309 
310 //VPunused Float_t vertexZ = mHiMicroEvent->VertexZ();
311 
312  for(Int_t i=0; i<nTrack; i++){
313  track =(StHiMicroTrack*) mHiMicroEvent->tracks()->At(i);
314 
315  //***************** spectra ***********************************
316 
317  if(!CutRc::AcceptTrackHalf(track)) {
318  continue; }
319  if(!CutRc::Accept(track)) continue;
320 
321  //**************************************************************
322  Float_t ptPr = track->PtPr();
323  Float_t eta = track->EtaPr();
324 //VPunused Float_t ptGl = track->PtGl();
325  //Float_t midZ = 100*TMath::Tan(track->DipAnglePr()) + vertexZ ;
326  //Float_t exitZ = 200*TMath::Tan(track->DipAnglePr()) + vertexZ ;
327  //
328  // efficiency correction here
329  //
330  Float_t eff = mEfficiencyMap->Eval(eta,ptPr);
331 
332  Float_t weight = (1./eff);
333 
334  //
335  // yield for both charge signs
336 
337  h2Yield->Fill(eta,ptPr);
338 
339  //
340  // counts
341  //
342  if(ptPr>1.5&&ptPr<6){
343  mNHiPtTrack++;
344  if(mNHiPtTrack%10000==0){
345  cout << "pt: " << ptPr << " eta: " << eta << " eff : " << eff << endl;
346  }
347  }
348  //
349  // spectra stuff here
350  //
351 
352  Int_t iCharge = (track->Charge()>0) ? 0 : 1; //plus is 0
353 
354  for(Int_t iBin=0; iBin<mNVarBin; iBin++){
355 
356  //
357  // both
358  //
359 
360  varBin[iBin].spec.h1Raw->Fill(ptPr);
361  varBin[iBin].mean.h1OneOverPt->Fill(ptPr,1./ptPr);
362  varBin[iBin].spec.h1EffCorrected->Fill(ptPr,weight);
363  varBin[iBin].spec.h1Corrected->Fill(ptPr,weight);
364 
365 
366  //
367  // plus/minus
368  //
369 
370  // mean stuff
371  //
372  varBin[iBin].meanPM[iCharge].h1OneOverPt->Fill(ptPr,1./ptPr);
373 
374  varBin[iBin].specPM[iCharge].h1Raw->Fill(ptPr);
375  varBin[iBin].specPM[iCharge].h1EffCorrected->Fill(ptPr,weight);
376  varBin[iBin].specPM[iCharge].h1Corrected->Fill(ptPr,weight);
377 
378  // vertex z east, west - make sure they stay on the same side
379  //
380 
381  } // varbin
382 
383  }
384 }
385 
386 //_____________________
387 
388 void
389 StHiSpectra::fillEventHistograms()
390 {
391  h1Centrality->Fill(mHiMicroEvent->Centrality());
392 
393  h1VertexZ->Fill(mHiMicroEvent->VertexZ());
394 
395  h1VertexZThin->Fill(mHiMicroEvent->VertexZ());
396 
397  if(CutRc::AcceptVertexZ(mHiMicroEvent))
398  h1CentralityCut->Fill(mHiMicroEvent->Centrality());
399 
400 // if(CutRc::AcceptZdcCtbCent(mHiMicroEvent)){
401 // h1VertexZCut->Fill(mHiMicroEvent->VertexZ());
402 // //h1VertexRCut->Fill(mHiMicroEvent->VertexR());
403 // }
404 }
405 
406 
407 //______________________
408 
409 Bool_t
410 StHiSpectra::acceptEvent(StHiMicroEvent* event)
411 {
412  return CutRc::Accept(event);
413 
414 }
415 
416 
417 //______________________
418 
419 void
420 StHiSpectra::finishHistograms()
421 {
422  h1NEvent->SetBinContent(1,mNEventAccepted);
423  h1EtaCut->SetBinContent(1,CutRc::mEta[0]);
424  h1EtaCut->SetBinContent(2,CutRc::mEta[1]);
425 
426  for(Int_t iBin=0; iBin<mNVarBin; iBin++){
427 
428  varBin[iBin].mean.h1WeightedMean->Divide(varBin[iBin].spec.h1Raw,varBin[iBin].mean.h1OneOverPt);
429 
430  for(Int_t iCharge=0; iCharge<2; iCharge++){
431 
432  varBin[iBin].meanPM[iCharge].h1WeightedMean->Divide(varBin[iBin].specPM[iCharge].h1Raw,varBin[iBin].meanPM[iCharge].h1OneOverPt);
433 
434  }
435  }
436 
437 }
438 
439 
440 
441 //_____________________________________
442 
443 ClassImp(StHiSpectra)
StHiSpectra::h1CentralityCut
TH1D * h1CentralityCut
just look at centrality stuff
Definition: StHiSpectra.h:65
StHiSpectra::h1EtaCut
TH1D * h1EtaCut
of events for the centrality class
Definition: StHiSpectra.h:74
Generated by   doxygen 1.8.5