StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructQAHists.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructQAHists.cxx,v 1.11 2012/11/16 21:19:07 prindle Exp $
4  *
5  * Author: Jeff Porter
6  *
7  **********************************************************************
8  *
9  * Description: base class for QA histogramming
10  *
11  ***********************************************************************/
12 #include "StEStructQAHists.h"
13 
14 #include "StEStructPool/EventMaker/StEStructEvent.h"
15 #include "StEStructPool/EventMaker/StEStructTrack.h"
16 #include "StEStructPool/EventMaker/StEStructCentrality.h"
17 #include "TH1F.h"
18 #include "TH1D.h"
19 #include "TH2F.h"
20 #include "TFile.h"
21 
22 ClassImp(StEStructQAHists)
23 
25  mEType=itype;
26  mCents[0]=mCents[1]=NULL;
27  mptCents[0]=mptCents[1]=mptCents[2]=NULL;
28  mRefMult=0;
29  mTotMult=mPosMult=mNegMult=0;
30  mTotMult4=mPosMult4=mNegMult4=0;
31  mntBins = 0;
32  mhasAIndex=false;
33  initHistograms();
34 };
35 
36 
37 StEStructQAHists::~StEStructQAHists(){ };
38 
39 //---------------------------------------------------------
40 //
41 // add your stuff below
42 //
43 //--------------------------------------------------------
44 void StEStructQAHists::initHistograms(){
45 
46  initBaseHistograms();
47  /* <-- add your analysis dependent set here --> */
48 
49 };
50 
51 void StEStructQAHists::fillHistograms(StEStructEvent* event, StEStructEventReader* reader){
52  fillBaseHistograms(event,reader);
53 };
54 
55 void StEStructQAHists::writeHistograms(TFile * tf){
56 
57  writeBaseHistograms(tf);
58  if(!mhasAIndex) writeTrackHistograms(tf);
59 
60  /* --- add your analysis dependent stuff here --- */
61 };
62 
63 
64 //
65 //
66 // ---- Below are common to all -----
67 //
68 //
69 
70 
71 void StEStructQAHists::initBaseHistograms(){
72 
73  StEStructCentrality* cent = StEStructCentrality::Instance();
74 
75  int mbNBins = cent->numCentralities();
76  if(mbNBins>0){
77  mCents[0] = new TH1D("cenClass","cenClass",mbNBins,-0.5,mbNBins-0.5);
78  mCents[1] = new TH1D("centralityDefs","centralityDefs",mbNBins,0.5,mbNBins+0.5);
79  }
80 
81  // fill the centrality mapping
82  for (int i=0;i<mbNBins;i++) mCents[1]->Fill(i,cent->centralityLimit(i));
83 
84  int mbNPtBins = cent->numPtCentralities();
85  if(mbNPtBins>0){
86  mptCents[0]= new TH1D("ptCenClass","ptCenClass",mbNPtBins,-0.5,mbNPtBins-0.5);
87  mptCents[1]= new TH1D("centralityPtDefs","centralityPtDefs",mbNPtBins,0.5,mbNPtBins+0.5);
88  }
89  for (int i=0;i<mbNPtBins;i++) mptCents[1]->Fill(i,cent->ptCentralityLimit(i));
90 
91 
92  int mbNPts = cent->numPts();
93  if(mbNPts>0){
94  mptCents[2] = new TH1D("ptRanges","ptRanges",mbNPts,-0.5,mbNPts+0.5);
95  }
96  for (int i=0;i<mbNPts;i++) mptCents[i]->Fill(i,cent->ptLimit(i));
97 
98  mRefMult = new TH1D("refMult","refMult for comparison",2500,-0.5,2500-0.5);
99  mTotMult = new TH1D("totalMultiplicity","totalMultiplicity",2500,-0.5,2500-0.5);
100  mPosMult = new TH1D("positiveMultiplicity","positiveMultiplicity",1500,-0.5,1500-0.5);
101  mNegMult = new TH1D("negativeMultiplicity","negativeMultiplicity",1500,-0.5,1500-0.5);
102  mTotMult4 = new TH1D("nTotOneQuarter","nTotOneQuarter",100,0.0,7.0);
103  mPosMult4 = new TH1D("nPosOneQuarter","nPosOneQuarter",100,0.0,7.0);
104  mNegMult4 = new TH1D("nNegOneQuarter","nNegOneQuarter",100,0.0,7.0);
105 
106  if(mEType==1) { // AA monte carlo generator
107 
108  aaGenImpact = new TH1D("impact","impact",100,0.0,20.0);
109  aaGen[0] = new TH2F("binary","binary",500,0.5,2000.5,50,0.0,20.0);
110  aaGen[1] = new TH2F("participant","participant",500,0.5,500.5,50,0.0,20.0);
111 
112  aaGenBin = new TH1D*[mbNBins-1];
113  aaGenPart = new TH1D*[mbNBins-1];
114  for(int i=0;i<mbNBins-1;i++){
115 
116  TString bName("binary_"); bName += i;
117  TString pName("participant_"); pName += i;
118  aaGenBin[i] = new TH1D(bName.Data(),bName.Data(),2000,0.5,2000.5);
119  aaGenPart[i] = new TH1D(pName.Data(),pName.Data(),2000,0.5,2000.5);
120 
121  }
122  }
123 
124  if(mEType==2){ //pp event generator ...
125 
126  ppELines = new TH1F*[mbNBins-1];
127  ppALines = new TH1F*[mbNBins-1];
128 
129  for(int i=0;i<mbNBins-1;i++){
130 
131  TString pName("partonLines_"); pName+=i;
132  ppELines[i] = new TH1F(pName.Data(),pName.Data(),1000,0.5,1000.5);
133  TString aName("allPartonLines_"); aName+=i;
134  ppALines[i] = new TH1F(aName.Data(),aName.Data(),200,0.,4000.);
135  }
136  }
137 
138 
139 };
140 
141 void StEStructQAHists::fillBaseHistograms(StEStructEvent* event, StEStructEventReader* reader){
142 
143  if(!event) return;
144 
145 
146  StEStructCentrality* cent = StEStructCentrality::Instance();
147  int ic = cent->centrality(event->Centrality());
148  mCents[0]->Fill(ic);
149  if(mptCents[0])mptCents[0]->Fill(ic);
150  mRefMult->Fill(event->RefMult());
151  mTotMult->Fill(event->Ntrack());
152  mPosMult->Fill(event->Npos());
153  mNegMult->Fill(event->Nneg());
154  mTotMult4->Fill(pow(event->Ntrack(),0.25));
155  mPosMult4->Fill(pow(event->Npos(),0.25));
156  mNegMult4->Fill(pow(event->Nneg(),0.25));
157 
158  if(mEType==1){
159  if(aaGenImpact)aaGenImpact->Fill(reader->getImpact());
160  if(aaGen[0])aaGen[0]->Fill(reader->getBinary(),reader->getImpact());
161  if(aaGen[1])aaGen[1]->Fill(reader->getParticipants(),reader->getImpact());
162  if(ic>=0){
163  if(aaGenBin[ic]) aaGenBin[ic]->Fill(reader->getBinary());
164  if(aaGenPart[ic])aaGenPart[ic]->Fill(reader->getParticipants());
165  }
166  } else if(mEType==2){
167  if(ic>=0){
168  ppELines[ic]->Fill(reader->getNPartonic());
169  ppALines[ic]->Fill(reader->getNPartonic());
170  }
171  }
172 
173 
174 };
175 
176 void StEStructQAHists::writeBaseHistograms(TFile* tf){
177  tf->cd();
178 
179  for(int i=0;i<2;i++)if(mCents[i])mCents[i]->Write();
180  for(int i=0;i<3;i++)if(mptCents[i])mptCents[i]->Write();
181  mRefMult->Write();
182  mTotMult->Write();
183  mPosMult->Write();
184  mNegMult->Write();
185  mTotMult4->Write();
186  mPosMult4->Write();
187  mNegMult4->Write();
188  if(mEType==1){
189  if(aaGenImpact)aaGenImpact->Write();
190  for(int i=0;i<2;i++)if(aaGen[i])aaGen[i]->Write();
191  for(int i=0;i<StEStructCentrality::Instance()->numCentralities()-1;i++){
192  if(aaGenBin[i])aaGenBin[i]->Write();
193  if(aaGenPart[i])aaGenPart[i]->Write();
194  }
195  } else if(mEType==2){
196  for(int i=0;i<StEStructCentrality::Instance()->numCentralities()-1;i++){
197  if(ppELines[i])ppELines[i]->Write();
198  if(ppALines[i])ppALines[i]->Write();
199  }
200  }
201 
202 };
203 
204 void StEStructQAHists::initTrackHistograms(int numBins, int aIndex){
205 
206  if(mntBins>0)return; // already been here in init
207 
208  mntBins = numBins;
209  int qbins = 2*numBins;
210 
211  mHEta = new TH1F*[qbins];
212  mHPhi = new TH1F*[qbins];
213  mHPt = new TH1F*[qbins];
214  mHYt = new TH1F*[qbins];
215  mHMass = new TH1F*[qbins];
216  mHdEdxPtot = new TH2F*[numBins];
217  mHToFPtot = new TH2F*[numBins];
218  mHEtaPt = new TH2F*[qbins];
219 
220  int nall = 100;
221  int nMass = 160;
222  float xetamin = -2.0;
223  float xetamax = 2.0;
224  float xphimin = -M_PI;
225  float xphimax = M_PI;
226  float xptmin = 0.;
227  float xptmax = 6.0;
228  float xytmin = 0.5;
229  float xytmax = 5.0;
230  float xmassmin = 0.0;
231  float xmassmax = 4.0;
232 
233  int nx = 101;
234  int ny = 101;
235  float xpmin = -2.5;
236  float xpmax = 2.5;
237  float ydmin = 0.;
238  float ydmax = 15.0e-06;
239  float ytofmin = 0.;
240  float ytofmax = 1.5;
241 
242  TString haIndex;
243  if(aIndex>=0){
244  mhasAIndex=true;
245  haIndex+="_A";
246  haIndex+=aIndex;
247  haIndex+="_";
248  }
249 
250  for(int i=0; i<numBins;i++){
251 
252  TString heta("Eta");
253  if(mhasAIndex)heta+=haIndex.Data();
254  if(numBins>1)heta+=i;
255  TString hpeta("Qp"); hpeta+=heta.Data();
256  mHEta[i]=new TH1F(hpeta.Data(),hpeta.Data(),nall,xetamin,xetamax);
257  TString hmeta("Qm"); hmeta+=heta.Data();
258  mHEta[i+numBins]=new TH1F(hmeta.Data(),hmeta.Data(),nall,xetamin,xetamax);
259 
260  TString hphi("Phi");
261  if(mhasAIndex)hphi+=haIndex.Data();
262  if(numBins>1)hphi+=i;
263  TString hpphi("Qp"); hpphi+=hphi.Data();
264  mHPhi[i]=new TH1F(hpphi.Data(),hpphi.Data(),nall,xphimin,xphimax);
265  TString hmphi("Qm"); hmphi+=hphi.Data();
266  mHPhi[i+numBins]=new TH1F(hmphi.Data(),hmphi.Data(),nall,xphimin,xphimax);
267 
268  TString hpt("Pt");
269  if(mhasAIndex)hpt+=haIndex.Data();
270  if(numBins>1)hpt+=i;
271  TString hppt("Qp"); hppt+=hpt.Data();
272  mHPt[i]=new TH1F(hppt.Data(),hppt.Data(),nall,xptmin,xptmax);
273  TString hmpt("Qm"); hmpt+=hpt.Data();
274  mHPt[i+numBins]=new TH1F(hmpt.Data(),hmpt.Data(),nall,xptmin,xptmax);
275 
276  TString hyt("Yt");
277  if(mhasAIndex)hyt+=haIndex.Data();
278  if(numBins>1)hyt+=i;
279  TString hpyt("Qp"); hpyt+=hyt.Data();
280  mHYt[i]=new TH1F(hpyt.Data(),hpyt.Data(),nall,xytmin,xytmax);
281  TString hmyt("Qm"); hmyt+=hyt.Data();
282  mHYt[i+numBins]=new TH1F(hmyt.Data(),hmyt.Data(),nall,xytmin,xytmax);
283 
284  TString hdedxP("dEdx_P");
285  if(mhasAIndex)hdedxP+=haIndex.Data();
286  if(numBins>1)hdedxP+=i;
287  mHdEdxPtot[i] = new TH2F(hdedxP.Data(),hdedxP.Data(),nx,xpmin,xpmax,ny,ydmin,ydmax);
288 
289  TString hmass("Mass");
290  if(mhasAIndex)hmass+=haIndex.Data();
291  if(numBins>1)hmass+=i;
292  TString hpmass("Qp"); hpmass+=hmass.Data();
293  mHMass[i]=new TH1F(hpmass.Data(),hpmass.Data(),nMass,xmassmin,xmassmax);
294  TString hmmass("Qm"); hmmass+=hmass.Data();
295  mHMass[i+numBins]=new TH1F(hmmass.Data(),hmmass.Data(),nMass,xmassmin,xmassmax);
296 
297  TString htofP("mass_P");
298  if(mhasAIndex)htofP+=haIndex.Data();
299  if(numBins>1)htofP+=i;
300  mHToFPtot[i] = new TH2F(htofP.Data(),htofP.Data(),nx,xpmin,xpmax,ny,ytofmin,ytofmax);
301 
302  TString hetapt("EtaPt");
303  if(mhasAIndex)hetapt+=haIndex.Data();
304  if(numBins>1)hetapt+=i;
305  TString hpetapt("Qp"); hpetapt+=hetapt.Data();
306  mHEtaPt[i]=new TH2F(hpetapt.Data(),hpetapt.Data(),20,-1.0,1.0,30,0.1,3.1);
307  TString hmetapt("Qm"); hmetapt+=hetapt.Data();
308  mHEtaPt[i+numBins]=new TH2F(hmetapt.Data(),hmetapt.Data(),20,-1.0,1.0,30,0.1,3.1);
309  }
310 
311 };
312 
313 
314 void StEStructQAHists::fillTrackHistograms(StEStructTrack* t, int ib) {
315 
316  if (mntBins==0 || ib>mntBins) {
317  return;
318  }
319  int i = ib;
320  if (t->Charge()<0) {
321  mHdEdxPtot[i]->Fill(-t->Ptot(),t->Dedx());
322  mHToFPtot[i]->Fill(-t->Ptot(),t->Mass());
323  i += mntBins;
324  } else {
325  mHdEdxPtot[i]->Fill(t->Ptot(),t->Dedx());
326  mHToFPtot[i]->Fill(t->Ptot(),t->Mass());
327  }
328 
329  mHEta[i]->Fill(t->Eta());
330  mHPhi[i]->Fill(t->Phi());
331  mHPt[i]->Fill(t->Pt());
332  mHYt[i]->Fill(t->Yt());
333  mHEtaPt[i]->Fill(t->Eta(),t->Pt());
334  mHMass[i]->Fill(t->Mass());
335 
336 };
337 
338 
339 void StEStructQAHists::writeTrackHistograms(TFile* tf){
340 
341  if(!tf || mntBins==0)return;
342  tf->cd();
343 
344  for(int i=0;i<2*mntBins;i++){
345  mHEta[i]->Write();
346  mHPhi[i]->Write();
347  mHPt[i]->Write();
348  mHYt[i]->Write();
349  mHEtaPt[i]->Write();
350  mHMass[i]->Write();
351  }
352  for(int i=0;i<mntBins;i++){
353  mHdEdxPtot[i]->Write();
354  mHToFPtot[i]->Write();
355  }
356 
357 };
358 
359 // Maybe implement these.
360 // Probably want to control via a doPairHistogram flag or something.
361 //void StEStructQAHists::initPairHistograms(int numBins, int aIndex){
362 //};
363 //void StEStructQAHists::fillPairHistograms(StEStructTrack* t1, StEStructTrack* t2, int ib, int after){
364 //};
365 //void StEStructQAHists::writePairHistograms(TFile* tf){
366 //};
367 /**********************************************************************
368  *
369  * $Log: StEStructQAHists.cxx,v $
370  * Revision 1.11 2012/11/16 21:19:07 prindle
371  * Moved EventCuts, TrackCuts to EventReader. Affects most readers.
372  * Added support to write and read EStructEvents.
373  * Cuts: 3D histo support, switch to control filling of histogram for reading EStructEvents
374  * EventCuts: A few new cuts
375  * MuDstReader: Add 2D to some histograms, treat ToFCut, PrimaryCuts, VertexRadius histograms like other cut histograms.
376  * QAHists: Add refMult
377  * TrackCuts: Add some hijing cuts.
378  *
379  * Revision 1.10 2011/08/02 20:31:25 prindle
380  * Change string handling
381  * Added event cuts for VPD, good fraction of global tracks are primary, vertex
382  * found only from tracks on single side of TPC, good fraction of primary tracks have TOF hits..
383  * Added methods to check if cuts imposed
384  * Added 2010 200GeV and 62 GeV, 2011 19 GeV AuAu datasets, 200 GeV pp2pp 2009 dataset.
385  * Added TOF vs. dEdx vs. p_t histograms
386  * Fix participant histograms in QAHists.
387  * Added TOFEMass cut in TrackCuts although I think we want to supersede this.
388  *
389  * Revision 1.9 2010/09/02 21:20:09 prindle
390  * Cuts: Add flag to not fill histograms. Important when scanning files for sorting.
391  * EventCuts: Add radius cut on vertex, ToF fraction cut. Merge 2004 AuAu 200 GeV datasets.
392  * Add 7, 11 and 39 GeV dataset selections
393  * MuDstReader: Add 2D histograms for vertex radius and ToF fraction cuts.
394  * Modify countGoodTracks to return number of dEdx and ToF pid identified tracks.
395  * Include code to set track pid information from Dst.
396  * QAHists: New ToF QA hists. Modify dEdx to include signed momentum.
397  *
398  * Revision 1.8 2010/06/23 22:29:47 prindle
399  * Hadd typo of 2004B instead of B2004 in EventCuts.cxx
400  * Added a couple of histograms in QAHists.
401  *
402  * Revision 1.7 2010/03/02 21:43:38 prindle
403  * Use outerHelix() for global tracks
404  * Add sensible triggerId histograms
405  * Starting to add support to sort events (available for Hijing)
406  *
407  * Revision 1.6 2008/03/19 22:02:00 prindle
408  * Updated some dataset definitions.
409  *
410  * Revision 1.5 2007/11/26 19:52:25 prindle
411  * Add cucu62, cucu200 2007ib production datasets.
412  * Included vertex cuts for case of ranked vertices. (Pass muDst pointer to EventCuts)
413  * Add n^(1/4) histograms to QAHists
414  *
415  * Revision 1.4 2006/04/27 22:20:10 prindle
416  * Some changes in trigger names for run periods.
417  * Changed a couple of the Hijing QA histograms.
418  *
419  * Revision 1.3 2006/04/11 17:50:50 prindle
420  * Remove inChain from constructor arguments (no longer used in macro)
421  *
422  * Revision 1.2 2006/04/06 00:54:03 prindle
423  * Tried to rationalize the way centrality is defined.
424  * Now the reader gives a float to StEStructEvent and this float is
425  * what is being used to define centrality. When we need a centrality
426  * bin index we pass this number into the centrality singleton object.
427  *
428  *
429  *
430  * Revision 1.1 2006/04/04 22:05:06 porter
431  * a handful of changes:
432  * - changed the StEStructAnalysisMaker to contain 1 reader not a list of readers
433  * - added StEStructQAHists object to contain histograms that did exist in macros or elsewhere
434  * - made centrality event cut taken from StEStructCentrality singleton
435  * - put in ability to get any max,min val from the cut class - one must call setRange in class
436  *
437  *
438  *********************************************************************/