StRoot: StFlowAnalysisMaker/StFlowScalarProdMaker.cxx Source File

00001 
00002 //
00003 // $Id: StFlowScalarProdMaker.cxx,v 1.12 2004/12/09 23:47:10 posk Exp $
00004 //
00005 // Authors: Method proposed by Art and Sergei, code written by Aihong
00006 //          Frame adopted from Art and Raimond's StFlowAnalysisMaker.
00008 //
00009 // Description:  Maker to analyze Flow using the scalar product method
00010 //
00012 
00013 #include <Stiostream.h>
00014 #include <stdlib.h>
00015 #include <math.h>
00016 #include <float.h>
00017 #include "StMaker.h"
00018 #include "StFlowScalarProdMaker.h"
00019 #include "StFlowMaker/StFlowMaker.h"
00020 #include "StFlowMaker/StFlowEvent.h"
00021 #include "StFlowMaker/StFlowConstants.h"
00022 #include "StFlowMaker/StFlowSelection.h"
00023 #include "StEnumerations.h"
00024 #include "PhysicalConstants.h"
00025 #include "SystemOfUnits.h"
00026 #include "TVector2.h"
00027 #include "TFile.h"
00028 #include "TString.h"
00029 #include "TH1.h"
00030 #include "TH2.h"
00031 #include "TProfile.h"
00032 #include "TProfile2D.h"
00033 #include "StMessMgr.h"
00034 #include "TMath.h"
00035 #define PR(x) cout << "##### FlowScalarProdAnalysis: " << (#x) << " = " << (x) << endl;
00036 
00037 ClassImp(StFlowScalarProdMaker)
00038 
00039 //-----------------------------------------------------------------------
00040 
00041 StFlowScalarProdMaker::StFlowScalarProdMaker(const Char_t* name): StMaker(name),
00042   MakerName(name) {
00043   pFlowSelect = new StFlowSelection();
00044   SetHistoRanges();
00045 }
00046 
00047 StFlowScalarProdMaker::StFlowScalarProdMaker(const Char_t* name,
00048     const StFlowSelection& flowSelect) : StMaker(name), MakerName(name) {
00049   pFlowSelect = new StFlowSelection(flowSelect); //copy constructor
00050   SetHistoRanges();
00051 }
00052 
00053 //-----------------------------------------------------------------------
00054 
00055 StFlowScalarProdMaker::~StFlowScalarProdMaker() {
00056 }
00057 
00058 //-----------------------------------------------------------------------
00059 
00060 Int_t StFlowScalarProdMaker::Make() {
00061   // Fill histograms
00062 
00063   // Get a pointer to StFlowEvent
00064   StFlowMaker* pFlowMaker = NULL;
00065   pFlowMaker = (StFlowMaker*)GetMaker("Flow");
00066   if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
00067   if (pFlowEvent) {
00068     if (pFlowSelect->Select(pFlowEvent)) {      // event selected
00069       FillFromFlowEvent();                      // get event quantities
00070       FillEventHistograms();                    // fill from FlowEvent
00071       if (pFlowEvent) FillParticleHistograms(); // fill particle histograms
00072     }
00073     if (Debug()) StMaker::PrintInfo();
00074   } else {
00075     gMessMgr->Info("##### FlowScalarProdMaker: FlowEvent pointer null");
00076   }
00077   
00078   return kStOK;
00079 }
00080 
00081 //-----------------------------------------------------------------------
00082 
00083 Int_t StFlowScalarProdMaker::Init() {
00084   // Book histograms
00085 
00086   float ptMaxPart = Flow::ptMaxPart;
00087   if (pFlowSelect->PtMaxPart()) {
00088     ptMaxPart = pFlowSelect->PtMaxPart();
00089   }
00090   int nPtBinsPart = Flow::nPtBinsPart;
00091   if (pFlowSelect->PtBinsPart()) {
00092     nPtBinsPart = pFlowSelect->PtBinsPart();
00093   }
00094   xLabel = "Pseudorapidity";
00095   if (strlen(pFlowSelect->PidPart()) != 0) { xLabel = "Rapidity"; }
00096 
00097   TString* histTitle;
00098 
00099   // for each selection
00100   for (int k = 0; k < Flow::nSels; k++) {
00101 
00102     // resolution
00103     histTitle = new TString("Flow_Res_ScalarProd_Sel");
00104     *histTitle += k+1;
00105     histFull[k].mHistRes = new TProfile(histTitle->Data(), histTitle->Data(),
00106       Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -1.*FLT_MAX, FLT_MAX, "");
00107     histFull[k].mHistRes->SetXTitle("Harmonic");
00108     histFull[k].mHistRes->SetYTitle("Resolution");
00109     delete histTitle;
00110 
00111     // vObs
00112     histTitle = new TString("Flow_vObs_ScalarProd_Sel");
00113     *histTitle += k+1;
00114     histFull[k].mHist_vObs = new TProfile(histTitle->Data(), histTitle->Data(),
00115       Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -10000., 10000., "");
00116     histFull[k].mHist_vObs->SetXTitle("Harmonic");
00117     histFull[k].mHist_vObs->SetYTitle("vObs (%)");
00118     delete histTitle;
00119     
00120     // for each harmonic
00121     for (int j = 0; j < Flow::nHars; j++) {
00122 
00123       // Flow observed
00124       histTitle = new TString("Flow_vObs2D_ScalarProd_Sel");
00125       *histTitle += k+1;
00126       histTitle->Append("_Har");
00127       *histTitle += j+1;
00128       histFull[k].histFullHar[j].mHist_vObs2D = new TProfile2D(histTitle->Data(),
00129         histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax, nPtBinsPart, 
00130                  Flow::ptMin, ptMaxPart, -10000., 10000., "");
00131       histFull[k].histFullHar[j].mHist_vObs2D->SetXTitle((char*)xLabel.Data());
00132       histFull[k].histFullHar[j].mHist_vObs2D->SetYTitle("Pt (GeV/c)");
00133       delete histTitle;
00134 
00135       // Flow observed profiles
00136       histTitle = new TString("Flow_vObsEta_ScalarProd_Sel");
00137       *histTitle += k+1;
00138       histTitle->Append("_Har");
00139       *histTitle += j+1;
00140       histFull[k].histFullHar[j].mHist_vObsEta = new TProfile(histTitle->Data(),
00141         histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax, 
00142                                                               -10000., 10000., "");
00143       histFull[k].histFullHar[j].mHist_vObsEta->SetXTitle((char*)xLabel.Data());
00144       histFull[k].histFullHar[j].mHist_vObsEta->SetYTitle("v (%)");
00145       delete histTitle;
00146 
00147       histTitle = new TString("Flow_vObsPt_ScalarProd_Sel");
00148       *histTitle += k+1;
00149       histTitle->Append("_Har");
00150       *histTitle += j+1;
00151       histFull[k].histFullHar[j].mHist_vObsPt = new TProfile(histTitle->Data(),
00152         histTitle->Data(), nPtBinsPart, Flow::ptMin, ptMaxPart, -10000., 10000., "");
00153       histFull[k].histFullHar[j].mHist_vObsPt->SetXTitle("Pt (GeV/c)");
00154       histFull[k].histFullHar[j].mHist_vObsPt->SetYTitle("v (%)");
00155       delete histTitle;
00156 
00157     }
00158   }
00159 
00160   gMessMgr->SetLimit("##### FlowScalarProd", 2);
00161   gMessMgr->Info("##### FlowScalarProdAnalysis: $Id: StFlowScalarProdMaker.cxx,v 1.12 2004/12/09 23:47:10 posk Exp $");
00162 
00163   return StMaker::Init();
00164 }
00165 
00166 
00167 
00168 //-----------------------------------------------------------------------
00169 
00170 void StFlowScalarProdMaker::FillFromFlowEvent() {
00171   // Get Q vectors from StFlowEvent
00172 
00173   for (int k = 0; k < Flow::nSels; k++) {
00174     pFlowSelect->SetSelection(k);
00175     for (int j = 0; j < Flow::nHars; j++) {
00176       pFlowSelect->SetHarmonic(j);
00177       for (int n = 0; n < Flow::nSubs; n++) {
00178         pFlowSelect->SetSubevent(n);
00179         int i = Flow::nSels*k + n;
00180         // sub-event quantities
00181         mQSub[i][j]=pFlowEvent->Q(pFlowSelect);
00182       }
00183       
00184       pFlowSelect->SetSubevent(-1);
00185       // full event quantities
00186       mQ[k][j]    = pFlowEvent->Q(pFlowSelect);
00187     }
00188   }
00189 
00190 }
00191 
00192 //-----------------------------------------------------------------------
00193 
00194 void StFlowScalarProdMaker::FillEventHistograms() {
00195   // The scaler product of the subevent Q vectors
00196 
00197   for (int k = 0; k < Flow::nSels; k++) {
00198     for (int j = 0; j < Flow::nHars; j++) {
00199       float order  = (float)(j+1);
00200 
00201       histFull[k].mHistRes->Fill(order, (mQSub[Flow::nSels*k + 0][j].X()) * 
00202           (mQSub[Flow::nSels*k + 1][j].X()) + (mQSub[Flow::nSels*k + 0][j].Y()) 
00203                                  * (mQSub[Flow::nSels*k + 1][j].Y()) ); 
00204 
00205     }
00206   }
00207 
00208 }
00209 
00210 //-----------------------------------------------------------------------
00211 
00212 void StFlowScalarProdMaker::FillParticleHistograms() {
00213   // Fill histograms from the particles
00214 
00215   // Initialize Iterator
00216   StFlowTrackCollection* pFlowTracks = pFlowEvent->TrackCollection();
00217   StFlowTrackIterator itr;
00218   
00219   for (itr = pFlowTracks->begin(); itr != pFlowTracks->end(); itr++) {
00220     StFlowTrack* pFlowTrack = *itr;
00221 
00222     float phi = pFlowTrack->Phi();
00223     if (phi < 0.) phi += twopi;
00224     float eta = pFlowTrack->Eta();
00225     float pt  = pFlowTrack->Pt();
00226     TVector2 q_i;
00227     TVector2 u_i;
00228 
00229     for (int k = 0; k < Flow::nSels; k++) {
00230       pFlowSelect->SetSelection(k);
00231       for (int j = 0; j < Flow::nHars; j++) {
00232         pFlowSelect->SetHarmonic(j);
00233 
00234         if (pFlowSelect->SelectPart(pFlowTrack)) {
00235           bool oddHar = (j+1) % 2;
00236           double order  = (double)(j+1);
00237           TVector2 mQ_i = mQ[k][j];
00238           
00239           // Remove autocorrelations
00240           if (pFlowSelect->Select(pFlowTrack)) {
00241             double phiWgt = pFlowEvent->PhiWeight(k, j, pFlowTrack);
00242             q_i.Set(phiWgt * cos(phi * order), phiWgt * sin(phi * order));
00243             mQ_i = mQ_i - q_i;
00244           }
00245                   
00246           // Caculate v for all particles selected
00247           u_i.Set(cos(phi*order), sin(phi*order));
00248           float v = (mQ_i.X()*u_i.X() + mQ_i.Y()*u_i.Y()) /perCent;
00249           float vFlip = v;
00250           if (eta < 0 && oddHar) vFlip *= -1;
00251           if (strlen(pFlowSelect->PidPart()) != 0) { 
00252             float rapidity = pFlowTrack->Y();
00253             histFull[k].histFullHar[j].mHist_vObs2D-> Fill(rapidity, pt, v);
00254             histFull[k].histFullHar[j].mHist_vObsEta->Fill(rapidity, v);
00255           } else {
00256             histFull[k].histFullHar[j].mHist_vObs2D-> Fill(eta, pt, v);
00257             histFull[k].histFullHar[j].mHist_vObsEta->Fill(eta, v);
00258           }
00259           histFull[k].histFullHar[j].mHist_vObsPt-> Fill(pt, vFlip);
00260           histFull[k].mHist_vObs->Fill(order, vFlip);
00261           
00262         }
00263       }
00264     }  
00265   }
00266 
00267 }
00268 
00269 
00270 //-----------------------------------------------------------------------
00271 
00272 Int_t StFlowScalarProdMaker::Finish() {
00273   // Calculates resolution and mean flow values
00274 
00275   TString* histTitle;
00276 
00277   double content;
00278   double error;
00279   double totalError;
00280 
00281   cout << endl << "##### Scalar Product Maker:" << endl;
00282 
00283   for (int k = 0; k < Flow::nSels; k++) {
00284 
00285     // Create the 1D v histogram
00286     histTitle = new TString("Flow_v_ScalarProd_Sel");
00287     *histTitle += k+1;
00288     histFull[k].mHist_v = 
00289       histFull[k].mHist_vObs->ProjectionX(histTitle->Data());
00290     histFull[k].mHist_v->SetTitle(histTitle->Data());
00291     histFull[k].mHist_v->SetXTitle("Harmonic");
00292     histFull[k].mHist_v->SetYTitle("v (%)");
00293     delete histTitle;
00294     AddHist(histFull[k].mHist_v);
00295 
00296     for (int j = 0; j < Flow::nHars; j++) {
00297 
00298       //Calculate the resolution
00299       mRes[k][j]    = ::sqrt(histFull[k].mHistRes->GetBinContent(j+1))*2.;
00300       mResErr[k][j] = (histFull[k].mHistRes->GetBinError(j+1))*2./mRes[k][j];
00301 
00302       // Create the v 2D histogram
00303       histTitle = new TString("Flow_v2D_ScalarProd_Sel");
00304       *histTitle += k+1;
00305       histTitle->Append("_Har");
00306       *histTitle += j+1;
00307       histFull[k].histFullHar[j].mHist_v2D = 
00308         histFull[k].histFullHar[j].mHist_vObs2D->ProjectionXY(histTitle->Data());
00309       histFull[k].histFullHar[j].mHist_v2D->SetTitle(histTitle->Data());
00310       histFull[k].histFullHar[j].mHist_v2D->SetXTitle((char*)xLabel.Data());
00311       histFull[k].histFullHar[j].mHist_v2D->SetYTitle("Pt (GeV/c)");
00312       histFull[k].histFullHar[j].mHist_v2D->SetZTitle("v (%)");
00313       delete histTitle;
00314       AddHist(histFull[k].histFullHar[j].mHist_v2D);
00315 
00316       // Create the 1D v histograms
00317       histTitle = new TString("Flow_vEta_ScalarProd_Sel");
00318       *histTitle += k+1;
00319       histTitle->Append("_Har");
00320       *histTitle += j+1;
00321       histFull[k].histFullHar[j].mHist_vEta = 
00322         histFull[k].histFullHar[j].mHist_vObsEta->ProjectionX(histTitle->Data());
00323       histFull[k].histFullHar[j].mHist_vEta->SetTitle(histTitle->Data());
00324       histFull[k].histFullHar[j].mHist_vEta->SetXTitle((char*)xLabel.Data());
00325       histFull[k].histFullHar[j].mHist_vEta->SetYTitle("v (%)");
00326       delete histTitle;
00327       AddHist(histFull[k].histFullHar[j].mHist_vEta);
00328 
00329       TString* histTitle = new TString("Flow_vPt_ScalarProd_Sel");
00330       *histTitle += k+1;
00331       histTitle->Append("_Har");
00332       *histTitle += j+1;
00333       histFull[k].histFullHar[j].mHist_vPt = 
00334         histFull[k].histFullHar[j].mHist_vObsPt->ProjectionX(histTitle->Data());
00335       histFull[k].histFullHar[j].mHist_vPt->SetTitle(histTitle->Data());
00336       histFull[k].histFullHar[j].mHist_vPt->SetXTitle("Pt (GeV/c)");
00337       histFull[k].histFullHar[j].mHist_vPt->SetYTitle("v (%)");
00338       delete histTitle;
00339       AddHist(histFull[k].histFullHar[j].mHist_vPt);
00340 
00341       // Calulate v = vObs / Resolution or Q.u/(2::sqrt(<Q_a.Q_b>))
00342       if (mRes[k][j]) {
00343         cout << "##### Resolution of the " << j+1 << "th harmonic = " << 
00344           mRes[k][j] << " +/- " << mResErr[k][j] << endl;
00345         // The systematic error of the resolution is not folded in.
00346         histFull[k].histFullHar[j].mHist_v2D ->Scale(1. / mRes[k][j]);
00347         histFull[k].histFullHar[j].mHist_vEta->Scale(1. / mRes[k][j]);
00348         histFull[k].histFullHar[j].mHist_vPt ->Scale(1. / mRes[k][j]);
00349         content = histFull[k].mHist_v->GetBinContent(j+1);
00350         content /=  mRes[k][j];
00351         histFull[k].mHist_v->SetBinContent(j+1, content);
00352         // The systematic error of the resolution is folded in.
00353         error = histFull[k].mHist_v->GetBinError(j+1);
00354         error /= mRes[k][j];
00355         totalError = fabs(content) * ::sqrt((error/content)*(error/content) +
00356                (mResErr[k][j]/mRes[k][j])*(mResErr[k][j]/mRes[k][j]));
00357         histFull[k].mHist_v->SetBinError(j+1, totalError);
00358         cout << "##### v" << j+1 << "= (" << content << " +/- " << error << 
00359           " +/- " << totalError << "(with syst.)) %" << endl;
00360       } else {
00361         cout << "##### Resolution of the " << j+1 << "th harmonic was zero."
00362              << endl;
00363         histFull[k].histFullHar[j].mHist_v2D ->Reset();
00364         histFull[k].histFullHar[j].mHist_vEta->Reset();
00365         histFull[k].histFullHar[j].mHist_vPt ->Reset();
00366         histFull[k].mHist_v->SetBinContent(j+1, 0.);
00367         histFull[k].mHist_v->SetBinError(j+1, 0.);
00368       }
00369     }
00370   }
00371   //GetHistList()->ls();
00372 
00373   // Write all histograms
00374   TFile histFile("flow.scalar.root", "RECREATE");
00375   GetHistList()->Write();
00376   histFile.Close();
00377   
00378   delete pFlowSelect;
00379 
00380   return StMaker::Finish();
00381 }
00382 
00383 //-----------------------------------------------------------------------
00384 
00385 void StFlowScalarProdMaker::SetHistoRanges(Bool_t ftpc_included) {
00386 
00387     if (ftpc_included) {
00388           mEtaMin = Flow::etaMin;
00389           mEtaMax = Flow::etaMax;
00390         mNEtaBins = Flow::nEtaBins;
00391     } else {
00392           mEtaMin = Flow::etaMinTpcOnly;
00393           mEtaMax = Flow::etaMaxTpcOnly;
00394         mNEtaBins = Flow::nEtaBinsTpcOnly;
00395     }
00396 
00397     return;
00398 }
00399 
00401 //
00402 // $Log: StFlowScalarProdMaker.cxx,v $
00403 // Revision 1.12  2004/12/09 23:47:10  posk
00404 // Minor changes in code formatting.
00405 // Added hist for TPC primary dca to AnalysisMaker.
00406 //
00407 // Revision 1.11  2003/09/02 17:58:11  perev
00408 // gcc 3.2 updates + WarnOff
00409 //
00410 // Revision 1.10  2003/07/07 21:58:20  posk
00411 // Made units of momentum GeV/c instead of GeV.
00412 //
00413 // Revision 1.9  2003/02/24 19:35:27  posk
00414 // Corrected mistake in the error of the resolution.
00415 // This only affected doubly integrated v values.
00416 //
00417 // Revision 1.8  2003/01/10 16:40:49  oldi
00418 // Several changes to comply with FTPC tracks:
00419 // - Switch to include/exclude FTPC tracks introduced.
00420 //   The same switch changes the range of the eta histograms.
00421 // - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
00422 // - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
00423 //   West, FarWest (depending on vertex.z()).
00424 // - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
00425 // - Cut to exclude mu-events with no primary vertex introduced.
00426 //   (This is possible for UPC events and FTPC tracks.)
00427 // - Global DCA cut for FTPC tracks added.
00428 // - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
00429 // - Charge cut for FTPC tracks added.
00430 //
00431 // Revision 1.7  2002/05/21 18:42:17  posk
00432 // Kirill's correction to minBias.C for bins with one count.
00433 //
00434 // Revision 1.6  2002/02/18 01:11:53  jeromel
00435 // Mandatory fix for FLT_MAX fix i.e. include float.h
00436 //
00437 // Revision 1.5  2002/02/13 22:31:46  posk
00438 // Pt Weight now also weights Phi Weight. Added Eta Weught, default=FALSE.
00439 //
00440 // Revision 1.4  2002/01/31 01:09:30  posk
00441 // *** empty log message ***
00442 //
00443 // Revision 1.3  2002/01/14 23:42:52  posk
00444 // Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
00445 //
00446 // Revision 1.2  2001/12/21 17:01:59  aihong
00447 // minor changes
00448 //
00449 // Revision 1.1  2001/12/18 23:46:47  aihong
00450 // install scalar product method
00451 //
00452 //