StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFlowScalarProdMaker.cxx
1 //
3 // $Id: StFlowScalarProdMaker.cxx,v 1.12 2004/12/09 23:47:10 posk Exp $
4 //
5 // Authors: Method proposed by Art and Sergei, code written by Aihong
6 // Frame adopted from Art and Raimond's StFlowAnalysisMaker.
8 //
9 // Description: Maker to analyze Flow using the scalar product method
10 //
12 
13 #include <Stiostream.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <float.h>
17 #include "StMaker.h"
18 #include "StFlowScalarProdMaker.h"
19 #include "StFlowMaker/StFlowMaker.h"
20 #include "StFlowMaker/StFlowEvent.h"
21 #include "StFlowMaker/StFlowConstants.h"
22 #include "StFlowMaker/StFlowSelection.h"
23 #include "StEnumerations.h"
24 #include "PhysicalConstants.h"
25 #include "SystemOfUnits.h"
26 #include "TVector2.h"
27 #include "TFile.h"
28 #include "TString.h"
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TProfile.h"
32 #include "TProfile2D.h"
33 #include "StMessMgr.h"
34 #include "TMath.h"
35 #define PR(x) cout << "##### FlowScalarProdAnalysis: " << (#x) << " = " << (x) << endl;
36 
37 ClassImp(StFlowScalarProdMaker)
38 
39 //-----------------------------------------------------------------------
40 
41 StFlowScalarProdMaker::StFlowScalarProdMaker(const Char_t* name): StMaker(name),
42  MakerName(name) {
43  pFlowSelect = new StFlowSelection();
44  SetHistoRanges();
45 }
46 
47 StFlowScalarProdMaker::StFlowScalarProdMaker(const Char_t* name,
48  const StFlowSelection& flowSelect) : StMaker(name), MakerName(name) {
49  pFlowSelect = new StFlowSelection(flowSelect); //copy constructor
50  SetHistoRanges();
51 }
52 
53 //-----------------------------------------------------------------------
54 
55 StFlowScalarProdMaker::~StFlowScalarProdMaker() {
56 }
57 
58 //-----------------------------------------------------------------------
59 
60 Int_t StFlowScalarProdMaker::Make() {
61  // Fill histograms
62 
63  // Get a pointer to StFlowEvent
64  StFlowMaker* pFlowMaker = NULL;
65  pFlowMaker = (StFlowMaker*)GetMaker("Flow");
66  if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
67  if (pFlowEvent) {
68  if (pFlowSelect->Select(pFlowEvent)) { // event selected
69  FillFromFlowEvent(); // get event quantities
70  FillEventHistograms(); // fill from FlowEvent
71  if (pFlowEvent) FillParticleHistograms(); // fill particle histograms
72  }
73  if (Debug()) StMaker::PrintInfo();
74  } else {
75  gMessMgr->Info("##### FlowScalarProdMaker: FlowEvent pointer null");
76  }
77 
78  return kStOK;
79 }
80 
81 //-----------------------------------------------------------------------
82 
83 Int_t StFlowScalarProdMaker::Init() {
84  // Book histograms
85 
86  float ptMaxPart = Flow::ptMaxPart;
87  if (pFlowSelect->PtMaxPart()) {
88  ptMaxPart = pFlowSelect->PtMaxPart();
89  }
90  int nPtBinsPart = Flow::nPtBinsPart;
91  if (pFlowSelect->PtBinsPart()) {
92  nPtBinsPart = pFlowSelect->PtBinsPart();
93  }
94  xLabel = "Pseudorapidity";
95  if (strlen(pFlowSelect->PidPart()) != 0) { xLabel = "Rapidity"; }
96 
97  TString* histTitle;
98 
99  // for each selection
100  for (int k = 0; k < Flow::nSels; k++) {
101 
102  // resolution
103  histTitle = new TString("Flow_Res_ScalarProd_Sel");
104  *histTitle += k+1;
105  histFull[k].mHistRes = new TProfile(histTitle->Data(), histTitle->Data(),
106  Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -1.*FLT_MAX, FLT_MAX, "");
107  histFull[k].mHistRes->SetXTitle("Harmonic");
108  histFull[k].mHistRes->SetYTitle("Resolution");
109  delete histTitle;
110 
111  // vObs
112  histTitle = new TString("Flow_vObs_ScalarProd_Sel");
113  *histTitle += k+1;
114  histFull[k].mHist_vObs = new TProfile(histTitle->Data(), histTitle->Data(),
115  Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -10000., 10000., "");
116  histFull[k].mHist_vObs->SetXTitle("Harmonic");
117  histFull[k].mHist_vObs->SetYTitle("vObs (%)");
118  delete histTitle;
119 
120  // for each harmonic
121  for (int j = 0; j < Flow::nHars; j++) {
122 
123  // Flow observed
124  histTitle = new TString("Flow_vObs2D_ScalarProd_Sel");
125  *histTitle += k+1;
126  histTitle->Append("_Har");
127  *histTitle += j+1;
128  histFull[k].histFullHar[j].mHist_vObs2D = new TProfile2D(histTitle->Data(),
129  histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax, nPtBinsPart,
130  Flow::ptMin, ptMaxPart, -10000., 10000., "");
131  histFull[k].histFullHar[j].mHist_vObs2D->SetXTitle((char*)xLabel.Data());
132  histFull[k].histFullHar[j].mHist_vObs2D->SetYTitle("Pt (GeV/c)");
133  delete histTitle;
134 
135  // Flow observed profiles
136  histTitle = new TString("Flow_vObsEta_ScalarProd_Sel");
137  *histTitle += k+1;
138  histTitle->Append("_Har");
139  *histTitle += j+1;
140  histFull[k].histFullHar[j].mHist_vObsEta = new TProfile(histTitle->Data(),
141  histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax,
142  -10000., 10000., "");
143  histFull[k].histFullHar[j].mHist_vObsEta->SetXTitle((char*)xLabel.Data());
144  histFull[k].histFullHar[j].mHist_vObsEta->SetYTitle("v (%)");
145  delete histTitle;
146 
147  histTitle = new TString("Flow_vObsPt_ScalarProd_Sel");
148  *histTitle += k+1;
149  histTitle->Append("_Har");
150  *histTitle += j+1;
151  histFull[k].histFullHar[j].mHist_vObsPt = new TProfile(histTitle->Data(),
152  histTitle->Data(), nPtBinsPart, Flow::ptMin, ptMaxPart, -10000., 10000., "");
153  histFull[k].histFullHar[j].mHist_vObsPt->SetXTitle("Pt (GeV/c)");
154  histFull[k].histFullHar[j].mHist_vObsPt->SetYTitle("v (%)");
155  delete histTitle;
156 
157  }
158  }
159 
160  gMessMgr->SetLimit("##### FlowScalarProd", 2);
161  gMessMgr->Info("##### FlowScalarProdAnalysis: $Id: StFlowScalarProdMaker.cxx,v 1.12 2004/12/09 23:47:10 posk Exp $");
162 
163  return StMaker::Init();
164 }
165 
166 
167 
168 //-----------------------------------------------------------------------
169 
170 void StFlowScalarProdMaker::FillFromFlowEvent() {
171  // Get Q vectors from StFlowEvent
172 
173  for (int k = 0; k < Flow::nSels; k++) {
174  pFlowSelect->SetSelection(k);
175  for (int j = 0; j < Flow::nHars; j++) {
176  pFlowSelect->SetHarmonic(j);
177  for (int n = 0; n < Flow::nSubs; n++) {
178  pFlowSelect->SetSubevent(n);
179  int i = Flow::nSels*k + n;
180  // sub-event quantities
181  mQSub[i][j]=pFlowEvent->Q(pFlowSelect);
182  }
183 
184  pFlowSelect->SetSubevent(-1);
185  // full event quantities
186  mQ[k][j] = pFlowEvent->Q(pFlowSelect);
187  }
188  }
189 
190 }
191 
192 //-----------------------------------------------------------------------
193 
194 void StFlowScalarProdMaker::FillEventHistograms() {
195  // The scaler product of the subevent Q vectors
196 
197  for (int k = 0; k < Flow::nSels; k++) {
198  for (int j = 0; j < Flow::nHars; j++) {
199  float order = (float)(j+1);
200 
201  histFull[k].mHistRes->Fill(order, (mQSub[Flow::nSels*k + 0][j].X()) *
202  (mQSub[Flow::nSels*k + 1][j].X()) + (mQSub[Flow::nSels*k + 0][j].Y())
203  * (mQSub[Flow::nSels*k + 1][j].Y()) );
204 
205  }
206  }
207 
208 }
209 
210 //-----------------------------------------------------------------------
211 
212 void StFlowScalarProdMaker::FillParticleHistograms() {
213  // Fill histograms from the particles
214 
215  // Initialize Iterator
216  StFlowTrackCollection* pFlowTracks = pFlowEvent->TrackCollection();
217  StFlowTrackIterator itr;
218 
219  for (itr = pFlowTracks->begin(); itr != pFlowTracks->end(); itr++) {
220  StFlowTrack* pFlowTrack = *itr;
221 
222  float phi = pFlowTrack->Phi();
223  if (phi < 0.) phi += twopi;
224  float eta = pFlowTrack->Eta();
225  float pt = pFlowTrack->Pt();
226  TVector2 q_i;
227  TVector2 u_i;
228 
229  for (int k = 0; k < Flow::nSels; k++) {
230  pFlowSelect->SetSelection(k);
231  for (int j = 0; j < Flow::nHars; j++) {
232  pFlowSelect->SetHarmonic(j);
233 
234  if (pFlowSelect->SelectPart(pFlowTrack)) {
235  bool oddHar = (j+1) % 2;
236  double order = (double)(j+1);
237  TVector2 mQ_i = mQ[k][j];
238 
239  // Remove autocorrelations
240  if (pFlowSelect->Select(pFlowTrack)) {
241  double phiWgt = pFlowEvent->PhiWeight(k, j, pFlowTrack);
242  q_i.Set(phiWgt * cos(phi * order), phiWgt * sin(phi * order));
243  mQ_i = mQ_i - q_i;
244  }
245 
246  // Caculate v for all particles selected
247  u_i.Set(cos(phi*order), sin(phi*order));
248  float v = (mQ_i.X()*u_i.X() + mQ_i.Y()*u_i.Y()) /perCent;
249  float vFlip = v;
250  if (eta < 0 && oddHar) vFlip *= -1;
251  if (strlen(pFlowSelect->PidPart()) != 0) {
252  float rapidity = pFlowTrack->Y();
253  histFull[k].histFullHar[j].mHist_vObs2D-> Fill(rapidity, pt, v);
254  histFull[k].histFullHar[j].mHist_vObsEta->Fill(rapidity, v);
255  } else {
256  histFull[k].histFullHar[j].mHist_vObs2D-> Fill(eta, pt, v);
257  histFull[k].histFullHar[j].mHist_vObsEta->Fill(eta, v);
258  }
259  histFull[k].histFullHar[j].mHist_vObsPt-> Fill(pt, vFlip);
260  histFull[k].mHist_vObs->Fill(order, vFlip);
261 
262  }
263  }
264  }
265  }
266 
267 }
268 
269 
270 //-----------------------------------------------------------------------
271 
272 Int_t StFlowScalarProdMaker::Finish() {
273  // Calculates resolution and mean flow values
274 
275  TString* histTitle;
276 
277  double content;
278  double error;
279  double totalError;
280 
281  cout << endl << "##### Scalar Product Maker:" << endl;
282 
283  for (int k = 0; k < Flow::nSels; k++) {
284 
285  // Create the 1D v histogram
286  histTitle = new TString("Flow_v_ScalarProd_Sel");
287  *histTitle += k+1;
288  histFull[k].mHist_v =
289  histFull[k].mHist_vObs->ProjectionX(histTitle->Data());
290  histFull[k].mHist_v->SetTitle(histTitle->Data());
291  histFull[k].mHist_v->SetXTitle("Harmonic");
292  histFull[k].mHist_v->SetYTitle("v (%)");
293  delete histTitle;
294  AddHist(histFull[k].mHist_v);
295 
296  for (int j = 0; j < Flow::nHars; j++) {
297 
298  //Calculate the resolution
299  mRes[k][j] = ::sqrt(histFull[k].mHistRes->GetBinContent(j+1))*2.;
300  mResErr[k][j] = (histFull[k].mHistRes->GetBinError(j+1))*2./mRes[k][j];
301 
302  // Create the v 2D histogram
303  histTitle = new TString("Flow_v2D_ScalarProd_Sel");
304  *histTitle += k+1;
305  histTitle->Append("_Har");
306  *histTitle += j+1;
307  histFull[k].histFullHar[j].mHist_v2D =
308  histFull[k].histFullHar[j].mHist_vObs2D->ProjectionXY(histTitle->Data());
309  histFull[k].histFullHar[j].mHist_v2D->SetTitle(histTitle->Data());
310  histFull[k].histFullHar[j].mHist_v2D->SetXTitle((char*)xLabel.Data());
311  histFull[k].histFullHar[j].mHist_v2D->SetYTitle("Pt (GeV/c)");
312  histFull[k].histFullHar[j].mHist_v2D->SetZTitle("v (%)");
313  delete histTitle;
314  AddHist(histFull[k].histFullHar[j].mHist_v2D);
315 
316  // Create the 1D v histograms
317  histTitle = new TString("Flow_vEta_ScalarProd_Sel");
318  *histTitle += k+1;
319  histTitle->Append("_Har");
320  *histTitle += j+1;
321  histFull[k].histFullHar[j].mHist_vEta =
322  histFull[k].histFullHar[j].mHist_vObsEta->ProjectionX(histTitle->Data());
323  histFull[k].histFullHar[j].mHist_vEta->SetTitle(histTitle->Data());
324  histFull[k].histFullHar[j].mHist_vEta->SetXTitle((char*)xLabel.Data());
325  histFull[k].histFullHar[j].mHist_vEta->SetYTitle("v (%)");
326  delete histTitle;
327  AddHist(histFull[k].histFullHar[j].mHist_vEta);
328 
329  TString* histTitle = new TString("Flow_vPt_ScalarProd_Sel");
330  *histTitle += k+1;
331  histTitle->Append("_Har");
332  *histTitle += j+1;
333  histFull[k].histFullHar[j].mHist_vPt =
334  histFull[k].histFullHar[j].mHist_vObsPt->ProjectionX(histTitle->Data());
335  histFull[k].histFullHar[j].mHist_vPt->SetTitle(histTitle->Data());
336  histFull[k].histFullHar[j].mHist_vPt->SetXTitle("Pt (GeV/c)");
337  histFull[k].histFullHar[j].mHist_vPt->SetYTitle("v (%)");
338  delete histTitle;
339  AddHist(histFull[k].histFullHar[j].mHist_vPt);
340 
341  // Calulate v = vObs / Resolution or Q.u/(2::sqrt(<Q_a.Q_b>))
342  if (mRes[k][j]) {
343  cout << "##### Resolution of the " << j+1 << "th harmonic = " <<
344  mRes[k][j] << " +/- " << mResErr[k][j] << endl;
345  // The systematic error of the resolution is not folded in.
346  histFull[k].histFullHar[j].mHist_v2D ->Scale(1. / mRes[k][j]);
347  histFull[k].histFullHar[j].mHist_vEta->Scale(1. / mRes[k][j]);
348  histFull[k].histFullHar[j].mHist_vPt ->Scale(1. / mRes[k][j]);
349  content = histFull[k].mHist_v->GetBinContent(j+1);
350  content /= mRes[k][j];
351  histFull[k].mHist_v->SetBinContent(j+1, content);
352  // The systematic error of the resolution is folded in.
353  error = histFull[k].mHist_v->GetBinError(j+1);
354  error /= mRes[k][j];
355  totalError = fabs(content) * ::sqrt((error/content)*(error/content) +
356  (mResErr[k][j]/mRes[k][j])*(mResErr[k][j]/mRes[k][j]));
357  histFull[k].mHist_v->SetBinError(j+1, totalError);
358  cout << "##### v" << j+1 << "= (" << content << " +/- " << error <<
359  " +/- " << totalError << "(with syst.)) %" << endl;
360  } else {
361  cout << "##### Resolution of the " << j+1 << "th harmonic was zero."
362  << endl;
363  histFull[k].histFullHar[j].mHist_v2D ->Reset();
364  histFull[k].histFullHar[j].mHist_vEta->Reset();
365  histFull[k].histFullHar[j].mHist_vPt ->Reset();
366  histFull[k].mHist_v->SetBinContent(j+1, 0.);
367  histFull[k].mHist_v->SetBinError(j+1, 0.);
368  }
369  }
370  }
371  //GetHistList()->ls();
372 
373  // Write all histograms
374  TFile histFile("flow.scalar.root", "RECREATE");
375  GetHistList()->Write();
376  histFile.Close();
377 
378  delete pFlowSelect;
379 
380  return StMaker::Finish();
381 }
382 
383 //-----------------------------------------------------------------------
384 
385 void StFlowScalarProdMaker::SetHistoRanges(Bool_t ftpc_included) {
386 
387  if (ftpc_included) {
388  mEtaMin = Flow::etaMin;
389  mEtaMax = Flow::etaMax;
390  mNEtaBins = Flow::nEtaBins;
391  } else {
392  mEtaMin = Flow::etaMinTpcOnly;
393  mEtaMax = Flow::etaMaxTpcOnly;
394  mNEtaBins = Flow::nEtaBinsTpcOnly;
395  }
396 
397  return;
398 }
399 
401 //
402 // $Log: StFlowScalarProdMaker.cxx,v $
403 // Revision 1.12 2004/12/09 23:47:10 posk
404 // Minor changes in code formatting.
405 // Added hist for TPC primary dca to AnalysisMaker.
406 //
407 // Revision 1.11 2003/09/02 17:58:11 perev
408 // gcc 3.2 updates + WarnOff
409 //
410 // Revision 1.10 2003/07/07 21:58:20 posk
411 // Made units of momentum GeV/c instead of GeV.
412 //
413 // Revision 1.9 2003/02/24 19:35:27 posk
414 // Corrected mistake in the error of the resolution.
415 // This only affected doubly integrated v values.
416 //
417 // Revision 1.8 2003/01/10 16:40:49 oldi
418 // Several changes to comply with FTPC tracks:
419 // - Switch to include/exclude FTPC tracks introduced.
420 // The same switch changes the range of the eta histograms.
421 // - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
422 // - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
423 // West, FarWest (depending on vertex.z()).
424 // - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
425 // - Cut to exclude mu-events with no primary vertex introduced.
426 // (This is possible for UPC events and FTPC tracks.)
427 // - Global DCA cut for FTPC tracks added.
428 // - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
429 // - Charge cut for FTPC tracks added.
430 //
431 // Revision 1.7 2002/05/21 18:42:17 posk
432 // Kirill's correction to minBias.C for bins with one count.
433 //
434 // Revision 1.6 2002/02/18 01:11:53 jeromel
435 // Mandatory fix for FLT_MAX fix i.e. include float.h
436 //
437 // Revision 1.5 2002/02/13 22:31:46 posk
438 // Pt Weight now also weights Phi Weight. Added Eta Weught, default=FALSE.
439 //
440 // Revision 1.4 2002/01/31 01:09:30 posk
441 // *** empty log message ***
442 //
443 // Revision 1.3 2002/01/14 23:42:52 posk
444 // Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
445 //
446 // Revision 1.2 2001/12/21 17:01:59 aihong
447 // minor changes
448 //
449 // Revision 1.1 2001/12/18 23:46:47 aihong
450 // install scalar product method
451 //
452 //
kStOK
Definition: Stypes.h:40
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
Generated by   doxygen 1.8.5