1 #include "StEpdEpFinder.h"
6 #include "TProfile2D.h"
7 #include "TClonesArray.h"
9 #include "StPicoEvent/StPicoEpdHit.h"
18 StEpdEpFinder::
StEpdEpFinder(
int nEventTypeBins,
char const* OutFileName,
char const* CorrectionFile) : mFormatUsed(2), mThresh(0.3), mMax(2.0), mWeightingScheme(0)
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";
28 mNumberOfEventTypeBins = nEventTypeBins;
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;
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;
51 for (
int ew=0; ew<2; ew++){
52 mPhiWeightInput[ew] = 0;
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));
61 mPhiWeightInput[ew] = (TH3D*)mCorrectionInputFile->Get(Form(
"PhiWeightEW%d",ew));
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));
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);
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);
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);
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);
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);
102 if (OrderOutsideRange(order)){
103 cout <<
"Order outside range - Ignoring your eta weights\n";
106 if (EtaWeights.GetYaxis()->GetNbins() != mNumberOfEventTypeBins){
107 cout <<
"Wrong number of EventType bins - Ignoring your eta weights\n";
110 mWeightingScheme = 1;
111 mEtaWeights[order-1] =
new TH2D;
112 EtaWeights.Copy(*mEtaWeights[order-1]);
116 if (OrderOutsideRange(order)){
117 cout <<
"Ignoring your ring weights\n";
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";
125 for (
int iring=0; iring<16; iring++){
126 mRingWeights[iring][order-1] = RingWeights[iring];
132 mCorrectionInputFile->Close();
134 for (
int ew=0; ew<2; ew++){
135 mPhiWeightOutput[ew]->Divide(mPhiAveraged[ew]);
136 delete mPhiAveraged[ew];
138 mCorrectionOutputFile->Write();
139 mCorrectionOutputFile->Close();
141 mResolutionOutputFile->Write();
142 mResolutionOutputFile->Close();
144 cout <<
"StEpdEpFinder out!\n\n";
150 if ((EventTypeId<0)||(EventTypeId>=mNumberOfEventTypeBins)){
151 cout <<
"You are asking for an undefined EventType - fail!\n";
157 double pi = TMath::Pi();
161 double TotalWeight4Ring[2][16][2];
162 double TotalWeight4Side[2][_EpOrderMax][2];
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;
168 for (
int iring=0; iring<16; iring++){
169 TotalWeight4Ring[iew][iring][phiWeightedOrNo] = 0;
176 for (
int hit=0;
hit<EpdHits->GetEntries();
hit++){
177 int tileId,ring,
TT,PP,EW,ADC;
179 switch(mFormatUsed) {
184 tileId = epdHit->
id();
186 ring = epdHit->
row();
190 nMip = epdHit->
nMIP();
195 std::cout <<
"You are requesting a format other than picoDst. It is not implemented\n";
196 std::cout <<
"YOU do it!\n";
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();
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);}
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);
220 double PhiWeightedTileWeight = TileWeight;
221 if (mPhiWeightInput[EW]) PhiWeightedTileWeight /= mPhiWeightInput[EW]->GetBinContent(PP,TT,EventTypeId+1);
222 TotalWeight4Ring[EW][ring-1][0] += TileWeight;
223 TotalWeight4Ring[EW][ring-1][1] += PhiWeightedTileWeight;
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;
228 TotalWeight4Side[EW][order-1][1] += fabs(etaWeight) * PhiWeightedTileWeight;
230 double Cosine = cos(phi*(
double)order);
231 double Sine = sin(phi*(
double)order);
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;
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;
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];
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];
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];
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];
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];
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];
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;
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);
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);
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);
327 for (
int order=1; order<_EpOrderMax+1; order++){
328 for (
int ewFull=0; ewFull<3; ewFull++){
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;
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;
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])));
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]));
368 double StEpdEpFinder::RingOrEtaWeight(
int ring,
double eta,
int order,
int EventTypeId){
369 switch (mWeightingScheme){
375 if (mEtaWeights[order-1]==0)
return 1.0;
376 int etaBin = mEtaWeights[order-1]->GetXaxis()->FindBin(fabs(eta));
377 return mEtaWeights[order-1]->GetBinContent(etaBin,EventTypeId+1);
381 return mRingWeights[ring-1][order-1];
389 double StEpdEpFinder::GetPsiInRange(
double Qx,
double Qy,
int order){
391 if ((Qx==0.0)||(Qy==0.0)) temp=-999;
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;
402 bool StEpdEpFinder::OrderOutsideRange(
int order){
404 cout <<
"\n*** Invalid order specified ***\n";
405 cout <<
"order must be 1 (for first-order plane) or higher\n";
408 if (order > _EpOrderMax) {
409 cout <<
"\n*** Invalid order specified ***\n";
410 cout <<
"Maximum order=" << _EpOrderMax <<
". To change, edit StEpdUtil/StEpdEpInfo.h\n";
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);
430 rep += Form(
"Weights for order %d W[j] means weight for EventType bin j\n",order);
432 for (
int EventTypeId=0; EventTypeId<mNumberOfEventTypeBins; EventTypeId++){rep += Form(
"\t\tW[%d]",EventTypeId);}
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));}
443 if (mWeightingScheme==2){
444 rep += Form(
"Here is the weighting: W[n] is weight for EP of order n\n");
446 for (
int order=1; order<=_EpOrderMax; order++){rep += Form(
"\tW[%d]",order);}
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]);}
Int_t tile() const
tile on the supersector [1,31]
void SetEtaWeights(int order, TH2D EtaWeight)
Int_t position() const
position of supersector on a wheel [1,12]
Int_t row() const
row number [1,16]
Int_t adc() const
ADC value [0,4095].
StEpdEpInfo Results(TClonesArray *EpdHits, TVector3 primVertex, int EventTypeID)
void SetRingWeights(int order, double *RingWeights)