StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEpdEpFinder.cxx
1 #include "StEpdEpFinder.h"
2 #include "TVector3.h"
3 #include "TH2D.h"
4 #include "TFile.h"
5 #include "TProfile.h"
6 #include "TProfile2D.h"
7 #include "TClonesArray.h"
8 #include "StEpdGeom.h"
9 #include "StPicoEvent/StPicoEpdHit.h"
10 #include "TMath.h"
11 #include <cassert>
12 #include <iostream>
13 //#include <fstream>
14 using namespace std;
15 
16 ClassImp(StEpdEpFinder)
17 
18 StEpdEpFinder::StEpdEpFinder(int nEventTypeBins, char const* OutFileName, char const* CorrectionFile) : mFormatUsed(2), mThresh(0.3), mMax(2.0), mWeightingScheme(0)
19 {
20 
21  cout << "\n**********\n* Welcome to the EPD Event Plane finder.\n"
22  << "* Please note that when you specify 'order' as an argument to a method,\n"
23  << "* then 1=first-order plane, 2=second-order plane, etc.\n"
24  << "* This code is currently configured to go up to order=" << _EpOrderMax << "\n"
25  << "* To include higher orders, edit StEpdUtil/StEpdEpInfo.h\n**********\n";
26 
27 
28  mNumberOfEventTypeBins = nEventTypeBins;
29 
30  mEpdGeom = new StEpdGeom();
31 
32  for (int iOrder=0; iOrder<_EpOrderMax; iOrder++){
33  mEtaWeights[iOrder] = 0;
34  for (int iRing=0; iRing<16; iRing++){
35  mRingWeights[iRing][iOrder]=1.0;
36  }
37  }
38 
39  // ----------------------------------- Stuff read from the "Correction File" -----------------------------------------
40  // "Shift correction" histograms that we INPUT and apply here
41  mCorrectionInputFile = new TFile(CorrectionFile,"READ");
42  if (mCorrectionInputFile->IsZombie()) {
43  std::cout << "Error opening file with Correction Histograms" << std::endl;
44  std::cout << "I will use no weighting at all." << std::endl;
45  for (int ewFull=0; ewFull<3; ewFull++){
46  for (int order=1; order<_EpOrderMax+1; order++){
47  mEpdShiftInput_sin[ewFull][order-1] = 0;
48  mEpdShiftInput_cos[ewFull][order-1] = 0;
49  }
50  }
51  for (int ew=0; ew<2; ew++){
52  mPhiWeightInput[ew] = 0;
53  }
54  }
55  else{
56  for (int ew=0; ew<2; ew++){
57  for (int order=1; order<_EpOrderMax+1; order++){
58  mEpdShiftInput_sin[ew][order-1] = (TProfile2D*)mCorrectionInputFile->Get(Form("EpdShiftEW%dPsi%d_sin",ew,order));
59  mEpdShiftInput_cos[ew][order-1] = (TProfile2D*)mCorrectionInputFile->Get(Form("EpdShiftEW%dPsi%d_cos",ew,order));
60  }
61  mPhiWeightInput[ew] = (TH3D*)mCorrectionInputFile->Get(Form("PhiWeightEW%d",ew));
62  }
63  for (int order=1; order<_EpOrderMax+1; order++){
64  mEpdShiftInput_sin[2][order-1] = (TProfile2D*)mCorrectionInputFile->Get(Form("EpdShiftFullEventPsi%d_sin",order));
65  mEpdShiftInput_cos[2][order-1] = (TProfile2D*)mCorrectionInputFile->Get(Form("EpdShiftFullEventPsi%d_cos",order));
66  }
67  }
68  // ----------------------------------- Stuff written to the "Correction File" -----------------------------------------
69  // "Shift correction" histograms that we produce and OUTPUT
70  mCorrectionOutputFile = new TFile(OutFileName,"RECREATE");
71  for (int ew=0; ew<2; ew++){
72  for (int order=1; order<_EpOrderMax+1; order++){
73  mEpdShiftOutput_sin[ew][order-1] = new TProfile2D(Form("EpdShiftEW%dPsi%d_sin",ew,order),Form("EpdShiftEW%dPsi%d_sin",ew,order),
74  _EpTermsMax,0.5,1.0*_EpTermsMax+.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
75  mEpdShiftOutput_cos[ew][order-1] = new TProfile2D(Form("EpdShiftEW%dPsi%d_cos",ew,order),Form("EpdShiftEW%dPsi%d_cos",ew,order),
76  _EpTermsMax,0.5,1.0*_EpTermsMax+.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
77  }
78  mPhiWeightOutput[ew] = new TH3D(Form("PhiWeightEW%d",ew),Form("Tile Weight divided by Averaged EW=%d",ew),12,0.5,12.5,31,0.5,31.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5); // bins are PP,TT,EventType
79  mPhiAveraged[ew] = new TH3D(Form("PhiAveraged%d",ew),Form("Average for this phi EW=%d",ew),12,0.5,12.5,31,0.5,31.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5); // just for normalization. discard after use
80  }
81  for (int order=1; order<_EpOrderMax+1; order++){
82  mEpdShiftOutput_sin[2][order-1] = new TProfile2D(Form("EpdShiftFullEventPsi%d_sin",order),
83  Form("EpdShiftFullEventPsi%d_sin",order),
84  _EpTermsMax,0.5,1.0*_EpTermsMax+.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
85  mEpdShiftOutput_cos[2][order-1] = new TProfile2D(Form("EpdShiftFullEventPsi%d_cos",order),
86  Form("EpdShiftFullEventPsi%d_cos",order),
87  _EpTermsMax,0.5,1.0*_EpTermsMax+.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
88  }
89 
90  // now stuff for the "Resolution File"
91  TString ResolutionFileName = "Resolution_";
92  ResolutionFileName += OutFileName;
93  mResolutionOutputFile = new TFile(ResolutionFileName,"RECREATE");
94  for (int order=1; order<=_EpOrderMax; order++){
95  mAveCosDeltaPsi[order-1] = new TProfile(Form("AveCosDeltaPsi%d",order),
96  Form("#LT cos(%d#Psi_{E,%d}-%d#Psi_{W,%d} #GT (weighted and shifted)",order,order,order,order),
97  nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
98  }
99 }
100 //------------------------------------------------------------------------------------------------------------
101 void StEpdEpFinder::SetEtaWeights(int order, TH2D EtaWeights){
102  if (OrderOutsideRange(order)){
103  cout << "Order outside range - Ignoring your eta weights\n";
104  return;
105  }
106  if (EtaWeights.GetYaxis()->GetNbins() != mNumberOfEventTypeBins){
107  cout << "Wrong number of EventType bins - Ignoring your eta weights\n";
108  return;
109  }
110  mWeightingScheme = 1;
111  mEtaWeights[order-1] = new TH2D;
112  EtaWeights.Copy(*mEtaWeights[order-1]);
113 }
114 
115 void StEpdEpFinder::SetRingWeights(int order, double* RingWeights){
116  if (OrderOutsideRange(order)){
117  cout << "Ignoring your ring weights\n";
118  return;
119  }
120  mWeightingScheme = 2;
121  if (order>_EpOrderMax){
122  std::cout << "You are setting weights for EP of order " << order << " which is too high! Increase _EpOrderMax in StEpdEpFinder.h\n";
123  assert(0);
124  }
125  for (int iring=0; iring<16; iring++){
126  mRingWeights[iring][order-1] = RingWeights[iring];
127  }
128 }
129 
131 
132  mCorrectionInputFile->Close();
133 
134  for (int ew=0; ew<2; ew++){
135  mPhiWeightOutput[ew]->Divide(mPhiAveraged[ew]);
136  delete mPhiAveraged[ew];
137  }
138  mCorrectionOutputFile->Write();
139  mCorrectionOutputFile->Close();
140 
141  mResolutionOutputFile->Write();
142  mResolutionOutputFile->Close();
143 
144  cout << "StEpdEpFinder out!\n\n";
145 }
146 
147 //==================================================================================================================
148 StEpdEpInfo StEpdEpFinder::Results(TClonesArray* EpdHits, TVector3 primVertex, int EventTypeId){
149 
150  if ((EventTypeId<0)||(EventTypeId>=mNumberOfEventTypeBins)){
151  cout << "You are asking for an undefined EventType - fail!\n";
152  assert(0);
153  }
154 
155  StEpdEpInfo result;
156 
157  double pi = TMath::Pi();
158 
159 
160  // This below is for normalizing Q-vectors
161  double TotalWeight4Ring[2][16][2]; // for normalizing Q-vector: indices EW, ring, (nonPhiWeighted or PhiWeighted)
162  double TotalWeight4Side[2][_EpOrderMax][2]; // for normalizing Q-vector: indices EW, order, (nonPhiWeighted or PhiWeighted) ** depends on Order because eta-weight depends on order
163  for (int iew=0; iew<2; iew++){
164  for (int phiWeightedOrNo=0; phiWeightedOrNo<2; phiWeightedOrNo++){
165  for (int order=1; order<_EpOrderMax+1; order++){
166  TotalWeight4Side[iew][order-1][phiWeightedOrNo] = 0;
167  }
168  for (int iring=0; iring<16; iring++){
169  TotalWeight4Ring[iew][iring][phiWeightedOrNo] = 0;
170  }
171  }
172  }
173 
174 
175  //--------------------------------- begin loop over hits ---------------------------------------
176  for (int hit=0; hit<EpdHits->GetEntries(); hit++){
177  int tileId,ring,TT,PP,EW,ADC;
178  float nMip;
179  switch(mFormatUsed) {
180  case 2 :
181  { // interesting. If a declaration is inside a case, it MUST be explicitly scoped with brackets {} or else error.
182  // https://en.cppreference.com/w/cpp/language/switch
183  StPicoEpdHit* epdHit = (StPicoEpdHit*)((*EpdHits)[hit]);
184  tileId = epdHit->id();
185  EW = (tileId<0)?0:1;
186  ring = epdHit->row();
187  TT = epdHit->tile();
188  PP = epdHit->position();
189  ADC = epdHit->adc();
190  nMip = epdHit->nMIP(); // malisa 20feb2019 - I have finally made the transition from ADC (next line) to truly nMip, now that calibrations are done.
191  // nMip = (TT<10)?(double)ADC/160.0:(double)ADC/115.0;
192  break;
193  }
194  default :
195  std::cout << "You are requesting a format other than picoDst. It is not implemented\n";
196  std::cout << "YOU do it!\n";
197  assert(0);
198  }
199  if (nMip<mThresh) continue;
200  double TileWeight = (nMip<mMax)?nMip:mMax;
201  TVector3 StraightLine = mEpdGeom->TileCenter(tileId) - primVertex;
202  double phi = StraightLine.Phi();
203  double eta = StraightLine.Eta();
204 
205  //---------------------------------
206  // fill Phi Weight histograms to be used in next iteration (if desired)
207  // Obviously, do this BEFORE phi weighting!
208  //---------------------------------
209  mPhiWeightOutput[EW]->Fill(PP,TT,EventTypeId,TileWeight);
210  if (TT==1){ for (int pp=1; pp<13; pp++) mPhiAveraged[EW]->Fill(pp,1,EventTypeId,TileWeight/12.0);}
211  else{
212  for (int pp=1; pp<13; pp++){
213  for (int tt=2*(ring-1); tt<2*ring; tt++) mPhiAveraged[EW]->Fill(pp,tt,EventTypeId,TileWeight/24.0);
214  }
215  }
216  //--------------------------------
217  // now calculate Q-vectors
218  //--------------------------------
219 
220  double PhiWeightedTileWeight = TileWeight;
221  if (mPhiWeightInput[EW]) PhiWeightedTileWeight /= mPhiWeightInput[EW]->GetBinContent(PP,TT,EventTypeId+1); // EventTypeId --> EventTypeId+1 on 27feb2020. Thanks to Xiaoyu Liu for finding this small bug. - MAL
222  TotalWeight4Ring[EW][ring-1][0] += TileWeight;
223  TotalWeight4Ring[EW][ring-1][1] += PhiWeightedTileWeight;
224 
225  for (int order=1; order<_EpOrderMax+1; order++){
226  double etaWeight = RingOrEtaWeight(ring,eta,order,EventTypeId);
227  TotalWeight4Side[EW][order-1][0] += fabs(etaWeight) * TileWeight; // yes the fabs() makes sense. The sign in the eta weight is equivalent to a trigonometric phase.
228  TotalWeight4Side[EW][order-1][1] += fabs(etaWeight) * PhiWeightedTileWeight; // yes the fabs() makes sense. The sign in the eta weight is equivalent to a trigonometric phase.
229 
230  double Cosine = cos(phi*(double)order);
231  double Sine = sin(phi*(double)order);
232 
233  result.QrawOneSide[EW][order-1][0] += etaWeight * TileWeight * Cosine;
234  result.QrawOneSide[EW][order-1][1] += etaWeight * TileWeight * Sine;
235  result.QringRaw[EW][order-1][0][ring-1] += TileWeight * Cosine;
236  result.QringRaw[EW][order-1][1][ring-1] += TileWeight * Sine;
237 
238  result.QphiWeightedOneSide[EW][order-1][0] += etaWeight * PhiWeightedTileWeight * Cosine;
239  result.QphiWeightedOneSide[EW][order-1][1] += etaWeight * PhiWeightedTileWeight * Sine;
240  result.QringPhiWeighted[EW][order-1][0][ring-1] += PhiWeightedTileWeight * Cosine;
241  result.QringPhiWeighted[EW][order-1][1][ring-1] += PhiWeightedTileWeight * Sine;
242  }
243  } // loop over hits
244  //--------------------------------- end loop over hits ---------------------------------------
245 
246  // Xiaoyu needs the weights used for each ring, so she can "un-normalize" the ring-by-ring Q-vectors. Okay...
247  for (int iew=0; iew<2; iew++){
248  for (int iring=1; iring<17; iring++){
249  result.RingSumWeightsRaw[iew][iring-1] = TotalWeight4Ring[iew][iring-1][0];
250  result.RingSumWeightsPhiWeighted[iew][iring-1] = TotalWeight4Ring[iew][iring-1][1];
251  }
252  for (int order=1; order<_EpOrderMax+1; order++){
253  result.WheelSumWeightsRaw[iew][order-1] = TotalWeight4Side[iew][order-1][0];
254  result.WheelSumWeightsPhiWeighted[iew][order-1] = TotalWeight4Side[iew][order-1][1];
255  }
256  }
257 
258 
259  // at this point, we are finished with the detector hits, and deal only with the Q-vectors,
260 
261  // first, normalize the Q's...
262  for (int EW=0; EW<2; EW++){
263  for (int order=1; order<_EpOrderMax+1; order++){
264  if (TotalWeight4Side[EW][order-1][0]>0.0001){
265  result.QrawOneSide[EW][order-1][0] /= TotalWeight4Side[EW][order-1][0];
266  result.QrawOneSide[EW][order-1][1] /= TotalWeight4Side[EW][order-1][0];
267  }
268  if (TotalWeight4Side[EW][order-1][1]>0.0001){
269  result.QphiWeightedOneSide[EW][order-1][0] /= TotalWeight4Side[EW][order-1][1];
270  result.QphiWeightedOneSide[EW][order-1][1] /= TotalWeight4Side[EW][order-1][1];
271  }
272  for (int ring=0; ring<16; ring++){
273  if (TotalWeight4Ring[EW][ring][0]>0.001){
274  result.QringRaw[EW][order-1][0][ring] /= TotalWeight4Ring[EW][ring][0];
275  result.QringRaw[EW][order-1][1][ring] /= TotalWeight4Ring[EW][ring][0];
276  }
277  if (TotalWeight4Ring[EW][ring][1]>0.001){
278  result.QringPhiWeighted[EW][order-1][0][ring] /= TotalWeight4Ring[EW][ring][1];
279  result.QringPhiWeighted[EW][order-1][1][ring] /= TotalWeight4Ring[EW][ring][1];
280  }
281  }
282  }
283  }
284 
285  // Before going any farther, flip the sign of the 1st-order Q-vector on the East side.
286  // I want the rapidity-odd first-order event plane. If you want something else, then
287  // write your own code, you weirdo ;-)
288  for (int xy=0; xy<2; xy++){
289  result.QrawOneSide[0][0][xy] *= -1.0;
290  result.QphiWeightedOneSide[0][0][xy] *= -1.0;
291  for (int iring=0; iring<16; iring++){
292  result.QringRaw[0][0][xy][iring] *= -1.0;
293  result.QringPhiWeighted[0][0][xy][iring] *= -1.0;
294  }
295  }
296 
297  // at this point, we are finished with the Q-vectors and just use them to get angles Psi
298 
299  //---------------------------------
300  // Calculate unshifted EP angles
301  // * note in this part, first we loop over East/West, then we do Full
302  //---------------------------------
303  for (int order=1; order<_EpOrderMax+1; order++){
304  for (int ew=0; ew<2; ew++){
305  result.PsiRaw[ew][order-1] = GetPsiInRange(result.QrawOneSide[ew][order-1][0],result.QrawOneSide[ew][order-1][1],order);
306  result.PsiPhiWeighted[ew][order-1] = GetPsiInRange(result.QphiWeightedOneSide[ew][order-1][0],result.QphiWeightedOneSide[ew][order-1][1],order);
307  for (int ring=1; ring<17; ring++){
308  result.PsiRingRaw[ew][order-1][ring-1] = GetPsiInRange(result.QringRaw[ew][order-1][0][ring-1],result.QringRaw[ew][order-1][1][ring-1],order);
309  result.PsiRingPhiWeighted[ew][order-1][ring-1] = GetPsiInRange(result.QringPhiWeighted[ew][order-1][0][ring-1],result.QringPhiWeighted[ew][order-1][1][ring-1],order);
310  } // loop over ring
311  } // loop over EW
312  // now the "Full Event" angles - these have the "ew index" = 2.
313  double Qx,Qy;
314  Qx = result.QrawOneSide[0][order-1][0] + result.QrawOneSide[1][order-1][0];
315  Qy = result.QrawOneSide[0][order-1][1] + result.QrawOneSide[1][order-1][1];
316  result.PsiRaw[2][order-1] = GetPsiInRange(Qx,Qy,order);
317 
318  Qx = result.QphiWeightedOneSide[0][order-1][0] + result.QphiWeightedOneSide[1][order-1][0];
319  Qy = result.QphiWeightedOneSide[0][order-1][1] + result.QphiWeightedOneSide[1][order-1][1];
320  result.PsiPhiWeighted[2][order-1] = GetPsiInRange(Qx,Qy,order);
321  } // loop over order
322 
323  //---------------------------------
324  // Now shift
325  //---------------------------------
326 
327  for (int order=1; order<_EpOrderMax+1; order++){
328  for (int ewFull=0; ewFull<3; ewFull++){ // 0,1,2 means East,West,Full
329  result.PsiPhiWeightedAndShifted[ewFull][order-1] = result.PsiPhiWeighted[ewFull][order-1];
330  if (mEpdShiftInput_sin[ewFull][order-1] != 0){
331  for (int i=1; i<=_EpTermsMax; i++){
332  double tmp = (double)(order*i);
333  double sinAve = mEpdShiftInput_sin[ewFull][order-1]->GetBinContent(i,EventTypeId+1);
334  double cosAve = mEpdShiftInput_cos[ewFull][order-1]->GetBinContent(i,EventTypeId+1);
335  result.PsiPhiWeightedAndShifted[ewFull][order-1] +=
336  2.0*(cosAve*sin(tmp*result.PsiPhiWeighted[ewFull][order-1]) - sinAve*cos(tmp*result.PsiPhiWeighted[ewFull][order-1]))/tmp;
337  }
338  double AngleWrapAround = 2.0*pi/(double)order;
339  if (result.PsiPhiWeightedAndShifted[ewFull][order-1]<0) result.PsiPhiWeightedAndShifted[ewFull][order-1] += AngleWrapAround;
340  else if (result.PsiPhiWeightedAndShifted[ewFull][order-1]>AngleWrapAround) result.PsiPhiWeightedAndShifted[ewFull][order-1] -= AngleWrapAround;
341  }
342  }
343  }
344 
345  // for the Resolutions file.
346  for (int order=1; order<=_EpOrderMax; order++){
347  mAveCosDeltaPsi[order-1]->Fill((double)EventTypeId,
348  cos(((double)order)*(result.PsiPhiWeightedAndShifted[0][order-1]-result.PsiPhiWeightedAndShifted[1][order-1])));
349  }
350 
351  //---------------------------------
352  // Now calculate shift histograms for a FUTURE run (if you want it)
353  //---------------------------------
354  for (int i=1; i<=_EpTermsMax; i++){
355  for (int ewFull=0; ewFull<3; ewFull++){
356  for (int order=1; order<_EpOrderMax+1; order++){
357  double tmp = (double)(order*i);
358  mEpdShiftOutput_sin[ewFull][order-1]->Fill(i,EventTypeId,sin(tmp*result.PsiPhiWeighted[ewFull][order-1]));
359  mEpdShiftOutput_cos[ewFull][order-1]->Fill(i,EventTypeId,cos(tmp*result.PsiPhiWeighted[ewFull][order-1]));
360  }
361  }
362  }
363 
364  return result;
365 }
366 
367 //==================================================================================================================
368 double StEpdEpFinder::RingOrEtaWeight(int ring, double eta, int order, int EventTypeId){
369  switch (mWeightingScheme){
370  case 0: // no weighting
371  return 1.0;
372  break;
373  case 1: // eta-based weights
374  { // must explicitly scope with {} if there is a declaration statement in a case
375  if (mEtaWeights[order-1]==0) return 1.0; // user never defined weights for this order
376  int etaBin = mEtaWeights[order-1]->GetXaxis()->FindBin(fabs(eta));
377  return mEtaWeights[order-1]->GetBinContent(etaBin,EventTypeId+1); // note the "+1" since EventTypeId starts at zero
378  break;
379  }
380  case 2: // ring-based weights
381  return mRingWeights[ring-1][order-1];
382  break;
383  default:
384  return 1.0;
385  }
386 }
387 
388 //==================================================================================================================
389 double StEpdEpFinder::GetPsiInRange(double Qx, double Qy, int order){
390  double temp;
391  if ((Qx==0.0)||(Qy==0.0)) temp=-999;
392  else{
393  temp = atan2(Qy,Qx)/((double)order);
394  double AngleWrapAround = 2.0*TMath::Pi()/(double)order;
395  if (temp<0.0) temp+= AngleWrapAround;
396  else if (temp>AngleWrapAround) temp -= AngleWrapAround;
397  }
398  return temp;
399 }
400 
401 //==================================================================================================================
402 bool StEpdEpFinder::OrderOutsideRange(int order){
403  if (order < 1) {
404  cout << "\n*** Invalid order specified ***\n";
405  cout << "order must be 1 (for first-order plane) or higher\n";
406  return true;
407  }
408  if (order > _EpOrderMax) {
409  cout << "\n*** Invalid order specified ***\n";
410  cout << "Maximum order=" << _EpOrderMax << ". To change, edit StEpdUtil/StEpdEpInfo.h\n";
411  return true;
412  }
413  return false;
414 }
415 
416 
417 //==================================================================================================================
419  TString rep = Form("This is the StEpdEpFinder Report:\n");
420  rep += Form("Number of EventType bins = %d\n",mNumberOfEventTypeBins);
421  rep += Form("Format of input data = %d [note: 0=StEvent, 1=StMuDst, 2=StPicoDst]\n",mFormatUsed);
422  rep += Form("Threshold (in MipMPV units) = %f and MAX weight = %f\n",mThresh,mMax);
423  rep += Form("Weighting scheme used=%d [note: 0=none, 1=eta-based weighting, 2=ring-based weighting]\n",mWeightingScheme);
424  if (mWeightingScheme==1){
425  for (int order=1; order<_EpOrderMax+1; order++){
426  if (mEtaWeights[order-1]==0){
427  rep += Form("No eta weighing for order n=%d\n",order);
428  }
429  else{
430  rep += Form("Weights for order %d W[j] means weight for EventType bin j\n",order);
431  rep += Form("eta");
432  for (int EventTypeId=0; EventTypeId<mNumberOfEventTypeBins; EventTypeId++){rep += Form("\t\tW[%d]",EventTypeId);}
433  rep += Form("\n");
434  for (int iEtaBin=1; iEtaBin<=mEtaWeights[order-1]->GetXaxis()->GetNbins(); iEtaBin++){
435  rep += Form("%f",mEtaWeights[order-1]->GetXaxis()->GetBinCenter(iEtaBin));
436  for (int EventTypeId=0; EventTypeId<mEtaWeights[order-1]->GetYaxis()->GetNbins(); EventTypeId++){rep += Form("\t%f",mEtaWeights[order-1]->GetBinContent(iEtaBin,EventTypeId+1));}
437  rep += Form("\n");
438  }
439  }
440  }
441  }
442  else{
443  if (mWeightingScheme==2){
444  rep += Form("Here is the weighting: W[n] is weight for EP of order n\n");
445  rep += Form("ring");
446  for (int order=1; order<=_EpOrderMax; order++){rep += Form("\tW[%d]",order);}
447  rep += Form("\n");
448  for (int iRing=1; iRing<=16; iRing++){
449  rep += Form("%d",iRing);
450  for (int order=1; order<=_EpOrderMax; order++){rep += Form("\t%f",mRingWeights[iRing-1][order-1]);}
451  rep += Form("\n");
452  }
453  }
454  }
455  return rep;
456 }
457 
Short_t id() const
Definition: StPicoEpdHit.h:90
Definition: T.h:18
Int_t tile() const
tile on the supersector [1,31]
Definition: StPicoEpdHit.h:94
TString Report()
void SetEtaWeights(int order, TH2D EtaWeight)
Int_t position() const
position of supersector on a wheel [1,12]
Definition: StPicoEpdHit.h:92
Int_t row() const
row number [1,16]
Definition: StPicoEpdHit.h:96
Int_t adc() const
ADC value [0,4095].
Definition: StPicoEpdHit.h:79
StEpdEpInfo Results(TClonesArray *EpdHits, TVector3 primVertex, int EventTypeID)
Float_t nMIP() const
Definition: StPicoEpdHit.h:104
void SetRingWeights(int order, double *RingWeights)