analysis/EpFinder/EpFinder.C

0001 #include "EpFinder.h"
0002 #include "TVector3.h"
0003 #include "TH2D.h"
0004 #include "TFile.h"
0005 #include "TProfile.h"
0006 #include "TProfile2D.h"
0007 #include "TClonesArray.h"
0008 #include "TMath.h"
0009 #include <cassert>
0010 #include <iostream>
0011 #include <vector>
0012
0013 using namespace std;
0014
0015
0016 EpFinder::EpFinder(int nEventTypeBins, char const* OutFileName, char const* CorrectionFile, int pbinsx, int pbinsy) : mThresh(0.3), mMax(2.0)
0017 {
0018
0019   cout << "\n**********\n*  Welcome to the Event Plane finder.\n"
0020        << "*  Please note that when you specify 'order' as an argument to a method,\n"
0021        << "*  then 1=first-order plane, 2=second-order plane, etc.\n"
0022        << "*  This code is currently configured to go up to order=" << _EpOrderMax << "\n";
0023   cout << "**********\n";
0024
0025   mNumberOfEventTypeBins = nEventTypeBins;
0026
0027   // ----------------------------------- Stuff read from the "Correction File" -----------------------------------------
0028   // "Shift correction" histograms that we INPUT and apply here
0029   mCorrectionInputFile = new TFile(CorrectionFile,"READ");
0030   if (mCorrectionInputFile->IsZombie()) {
0031     std::cout << "EPFinder: Error opening file with Correction Histograms" << std::endl;
0032     std::cout << "EPFinder: I will use no weighting at all." << std::endl;
0033     for (int order=1; order<_EpOrderMax+1; order++){
0034       mEpShiftInput_sin[order-1] = NULL;
0035       mEpShiftInput_cos[order-1] = NULL;
0036     }
0037     mPhiWeightInput = NULL;
0038   }
0039   else{
0040     for (int order=1; order<_EpOrderMax+1; order++){
0041       mEpShiftInput_sin[order-1] = (TProfile2D*)mCorrectionInputFile->Get(Form("EpShiftPsi%d_sin",order));
0042       mEpShiftInput_cos[order-1] = (TProfile2D*)mCorrectionInputFile->Get(Form("EpShiftPsi%d_cos",order));
0043     }
0044     mPhiWeightInput = (TH2D*)mCorrectionInputFile->Get(Form("PhiWeight"));
0045   }
0046   // ----------------------------------- Stuff written to the "Correction File" -----------------------------------------
0047   // "Shift correction" histograms that we produce and OUTPUT
0048   mCorrectionOutputFile = new TFile(OutFileName,"RECREATE");
0049   for (int order=1; order<_EpOrderMax+1; order++){
0050     mEpShiftOutput_sin[order-1] = new TProfile2D(Form("EpShiftPsi%d_sin",order),Form("EpShiftPsi%d_sin",order),
0051                           _EpTermsMax,0.5,1.0*_EpTermsMax+.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
0052     mEpShiftOutput_cos[order-1] = new TProfile2D(Form("EpShiftPsi%d_cos",order),Form("EpShiftPsi%d_cos",order),
0053                           _EpTermsMax,0.5,1.0*_EpTermsMax+.5,nEventTypeBins,-0.5,(double)nEventTypeBins-0.5,-1.0,1.0);
0054   }
0055   // Phi weighting histograms based on requested binning
0056   if(pbinsx<=0) pbinsx = 1;
0057   if(pbinsy<=0) pbinsy = 1;
0058
0059   // binning tuned for FEMC, to be made generic in future JGL 8/28/2019
0060   // bins are ix,iy
0061   mPhiWeightOutput   = new TH2D(Form("PhiWeight"),Form("Phi Weight"),pbinsx,-0.5,(pbinsx-0.5),pbinsy,-0.5,(pbinsy-0.5));
0062   // just for normalization. discard after use
0063   mPhiAveraged       = new TH2D(Form("PhiAveraged"),Form("Phi Average"),pbinsx,-0.5,(pbinsx-0.5),pbinsy,-0.5,(pbinsy-0.5));
0064
0065 }
0066
0067 void EpFinder::Finish(){
0068
0069   mCorrectionInputFile->Close();
0070
0071   mPhiWeightOutput->Divide(mPhiAveraged);
0072   delete mPhiAveraged;
0073
0074   mCorrectionOutputFile->Write();
0075   mCorrectionOutputFile->Close();
0076
0077   cout << "EpFinder is finished!\n\n";
0078 }
0079
0080 //==================================================================================================================
0081 EpInfo EpFinder::Results(std::vector<EpHit> *EpdHits, int EventTypeId){
0082
0083   if ((EventTypeId<0)||(EventTypeId>=mNumberOfEventTypeBins)){
0084     cout << "You are asking for an undefined EventType - fail!\n";
0085     assert(0);
0086   }
0087
0088   EpInfo result;
0089
0090   double pi = TMath::Pi();
0091
0092   // This below is for normalizing Q-vectors
0093   double TotalWeight4Side[_EpOrderMax][2];       // for normalizing Q-vector: order, (nonPhiWeighted or PhiWeighted)  ** depends on Order because eta-weight depends on order
0094   for (int phiWeightedOrNo=0; phiWeightedOrNo<2; phiWeightedOrNo++){
0095     for (int order=1; order<_EpOrderMax+1; order++){
0096       TotalWeight4Side[order-1][phiWeightedOrNo] = 0;
0097     }
0098   }
0099
0100   //--------------------------------- begin loop over hits ---------------------------------------
0101   for (unsigned int hit=0; hit<EpdHits->size(); hit++){
0102
0103     float nMip = EpdHits->at(hit).nMip;
0104     if (nMip<mThresh) continue;
0105     double TileWeight = (nMip<mMax)?nMip:mMax;
0106     int idx_x =  EpdHits->at(hit).ix;
0107     int idx_y =  EpdHits->at(hit).iy;
0108     double phi = EpdHits->at(hit).phi;
0109
0110     //---------------------------------
0111     // fill Phi Weight histograms to be used in next iteration (if desired)
0112     // Obviously, do this BEFORE phi weighting!
0113     //---------------------------------
0114
0115     if(EpdHits->at(hit).samePhi){
0116       mPhiWeightOutput->Fill(idx_x,idx_y,TileWeight);
0117       for(unsigned int i = 0; i<EpdHits->at(hit).samePhi->size(); i++){
0118     float x = EpdHits->at(hit).samePhi->at(i).first;
0119     float y = EpdHits->at(hit).samePhi->at(i).second;
0120     mPhiAveraged->Fill(x,y,TileWeight/EpdHits->at(hit).samePhi->size());
0121       }
0122     }
0123
0124     //--------------------------------
0125     // now calculate Q-vectors
0126     //--------------------------------
0127
0128     double PhiWeightedTileWeight = TileWeight;
0129     if (mPhiWeightInput) PhiWeightedTileWeight /= mPhiWeightInput->GetBinContent(idx_x+1,idx_y+1);
0130
0131     for (int order=1; order<_EpOrderMax+1; order++){
0132       double etaWeight = 1.0; // not implemented - JGL 8/27/2019
0133       TotalWeight4Side[order-1][0] += fabs(etaWeight) * TileWeight;             // yes the fabs() makes sense.  The sign in the eta weight is equivalent to a trigonometric phase.
0134       TotalWeight4Side[order-1][1] += fabs(etaWeight) * PhiWeightedTileWeight;  // yes the fabs() makes sense.  The sign in the eta weight is equivalent to a trigonometric phase.
0135
0136       double Cosine = cos(phi*(double)order);
0137       double Sine   = sin(phi*(double)order);
0138
0139       result.QrawOneSide[order-1][0]      += etaWeight * TileWeight * Cosine;
0140       result.QrawOneSide[order-1][1]      += etaWeight * TileWeight * Sine;
0141
0142       result.QphiWeightedOneSide[order-1][0]      += etaWeight * PhiWeightedTileWeight * Cosine;
0143       result.QphiWeightedOneSide[order-1][1]      += etaWeight * PhiWeightedTileWeight * Sine;
0144
0145     }
0146   }  // loop over hits
0147
0148   //--------------------------------- end loop over hits ---------------------------------------
0149
0150   // Weights used, so you can "un-normalize" the ring-by-ring Q-vectors.
0151   for (int order=1; order<_EpOrderMax+1; order++){
0152     result.WheelSumWeightsRaw[order-1]         = TotalWeight4Side[order-1][0];
0153     result.WheelSumWeightsPhiWeighted[order-1] = TotalWeight4Side[order-1][1];
0154   }
0155
0156   // at this point, we are finished with the detector hits, and deal only with the Q-vectors,
0157
0158   // first, normalize the Q's...
0159   for (int order=1; order<_EpOrderMax+1; order++){
0160     if (TotalWeight4Side[order-1][0]>0.0001){
0161       result.QrawOneSide[order-1][0] /= TotalWeight4Side[order-1][0];
0162       result.QrawOneSide[order-1][1] /= TotalWeight4Side[order-1][0];
0163     }
0164     if (TotalWeight4Side[order-1][1]>0.0001){
0165       result.QphiWeightedOneSide[order-1][0] /= TotalWeight4Side[order-1][1];
0166       result.QphiWeightedOneSide[order-1][1] /= TotalWeight4Side[order-1][1];
0167     }
0168   }
0169
0170   // at this point, we are finished with the Q-vectors and just use them to get angles Psi
0171
0172   //---------------------------------
0173   // Calculate unshifted EP angles
0174   //---------------------------------
0175   for (int order=1; order<_EpOrderMax+1; order++){
0176     result.PsiRaw[order-1]                       = GetPsiInRange(result.QrawOneSide[order-1][0],result.QrawOneSide[order-1][1],order);
0177     result.PsiPhiWeighted[order-1]               = GetPsiInRange(result.QphiWeightedOneSide[order-1][0],result.QphiWeightedOneSide[order-1][1],order);
0178   } // loop over order
0179
0180   //---------------------------------
0181   // Now shift
0182   //---------------------------------
0183   for (int order=1; order<_EpOrderMax+1; order++){
0184     result.PsiPhiWeightedAndShifted[order-1] = result.PsiPhiWeighted[order-1];
0185     if (mEpShiftInput_sin[order-1]) {
0186       for (int i=1; i<=_EpTermsMax; i++){
0187     double tmp = (double)(order*i);
0188     double sinAve = mEpShiftInput_sin[order-1]->GetBinContent(i,EventTypeId+1);    /// note the "+1" since EventTypeId begins at zero
0189     double cosAve = mEpShiftInput_cos[order-1]->GetBinContent(i,EventTypeId+1);    /// note the "+1" since EventTypeId begins at zero
0190     result.PsiPhiWeightedAndShifted[order-1] +=
0191       2.0*(cosAve*sin(tmp*result.PsiPhiWeighted[order-1]) - sinAve*cos(tmp*result.PsiPhiWeighted[order-1]))/tmp;
0192       }
0193       double AngleWrapAround = 2.0*pi/(double)order;
0194       if (result.PsiPhiWeightedAndShifted[order-1]<0) result.PsiPhiWeightedAndShifted[order-1] += AngleWrapAround;
0195       else if (result.PsiPhiWeightedAndShifted[order-1]>AngleWrapAround) result.PsiPhiWeightedAndShifted[order-1] -= AngleWrapAround;
0196     }
0197   }
0198
0199   //---------------------------------
0200   // Now calculate shift histograms for a FUTURE run (if you want it)
0201   //---------------------------------
0202   for (int i=1; i<=_EpTermsMax; i++){
0203     for (int order=1; order<_EpOrderMax+1; order++){
0204       double tmp = (double)(order*i);
0205       mEpShiftOutput_sin[order-1]->Fill(i,EventTypeId,sin(tmp*result.PsiPhiWeighted[order-1]));
0206       mEpShiftOutput_cos[order-1]->Fill(i,EventTypeId,cos(tmp*result.PsiPhiWeighted[order-1]));
0207     }
0208   }
0209
0210   return result;
0211 }
0212
0213 double EpFinder::GetPsiInRange(double Qx, double Qy, int order){
0214   double temp;
0215   if ((Qx==0.0)||(Qy==0.0)) temp=-999;
0216   else{
0217     temp = atan2(Qy,Qx)/((double)order);
0218     double AngleWrapAround = 2.0*TMath::Pi()/(double)order;
0219     if (temp<0.0) temp+= AngleWrapAround;
0220     else if (temp>AngleWrapAround) temp -= AngleWrapAround;
0221   }
0222   return temp;
0223 }
0224
0225 bool EpFinder::OrderOutsideRange(int order){
0226   if (order < 1) {
0227     cout << "\n*** Invalid order specified ***\n";
0228     cout << "order must be 1 (for first-order plane) or higher\n";
0229     return true;
0230   }
0231   if (order > _EpOrderMax) {
0232     cout << "\n*** Invalid order specified ***\n";
0233     cout << "Maximum order=" << _EpOrderMax << ". To change, edit StEpdUtil/StEpdEpInfo.h\n";
0234     return true;
0235   }
0236   return false;
0237 }
0238
0239 TString EpFinder::Report(){
0240   TString rep = Form("This is the EpFinder Report:\n");
0241   rep += Form("Number of EventType bins = %d\n",mNumberOfEventTypeBins);
0242   rep += Form("Threshold (in MipMPV units) = %f  and MAX weight = %f\n",mThresh,mMax);
0243   return rep;
0244 }
0245