StRoot: StFlowAnalysisMaker/StFlowCumulantMaker.cxx Source File

00001 
00002 //
00003 // $Id: StFlowCumulantMaker.cxx,v 1.22 2006/02/22 19:12:29 posk Exp $
00004 //
00005 // Authors:  Aihong Tang, Kent State U. Oct 2001
00006 //           Frame adopted from Art and Raimond's StFlowAnalysisMaker.
00007 //
00009 //
00010 // Description:  Maker to analyze Flow using the latest cumulant method.
00011 //                refer to Phy. Rev. C63 (2001) 054906 (old new method)
00012 //                and      Phy. Rev. C64 (2001) 054901 (new new method)
00013 //                and      nucl-ex/0110016             (Practical Guide)
00014 //                and      PRy. Rev. C66 (2002) 014905 (directed flow from elliptic flow v1{3})
00015 //                Eq. numbers are from the PG or new new if not specified.
00016 //
00017 //               Anything that has "Mix" in it is for v1{3} calculation 
00018 //
00020 
00021 #include <Stiostream.h>
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #include <float.h>
00025 #include "StMaker.h"
00026 #include "StFlowCumulantMaker.h"
00027 #include "StFlowMaker/StFlowMaker.h"
00028 #include "StFlowMaker/StFlowEvent.h"
00029 #include "StFlowMaker/StFlowConstants.h"
00030 #include "StFlowMaker/StFlowSelection.h"
00031 #include "StEnumerations.h"
00032 #include "PhysicalConstants.h"
00033 #include "SystemOfUnits.h"
00034 #include "TFile.h"
00035 #include "TString.h"
00036 #include "TObjString.h"
00037 #include "TVectorD.h"
00038 #include "TH1.h"
00039 #include "TH2.h"
00040 #include "TProfile.h"
00041 #include "TProfile2D.h"
00042 #include "TOrdCollection.h"
00043 #include "StMessMgr.h"
00044 #include "TMath.h"
00045 
00046 #define PR(x) cout << "##### FlowCumulantAnalysis: " << (#x) << " = " << (x) << endl;
00047 
00048 ClassImp(StFlowCumulantMaker)
00049 
00050 //-----------------------------------------------------------------------
00051 
00052 StFlowCumulantMaker::StFlowCumulantMaker(const Char_t* name): StMaker(name),
00053   MakerName(name) {
00054   pFlowSelect = new StFlowSelection();
00055   SetHistoRanges();
00056 }
00057 
00058 StFlowCumulantMaker::StFlowCumulantMaker(const Char_t* name,
00059                                          const StFlowSelection& flowSelect) :
00060   StMaker(name), MakerName(name) {
00061   pFlowSelect = new StFlowSelection(flowSelect); // copy constructor
00062   SetHistoRanges();
00063 }
00064 
00065 //-----------------------------------------------------------------------
00066 
00067 StFlowCumulantMaker::~StFlowCumulantMaker() {
00068 }
00069 
00070 //-----------------------------------------------------------------------
00071 
00072 Int_t StFlowCumulantMaker::Make() {
00073   // Fill histograms
00074 
00075   // Get a pointer to StFlowEvent
00076   StFlowMaker* pFlowMaker = NULL;
00077   pFlowMaker = (StFlowMaker*)GetMaker("Flow");
00078   if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
00079   if (pFlowEvent && pFlowSelect->Select(pFlowEvent)) { // event selected
00080     FillFromFlowEvent();                               // get event quantities
00081     FillEventHistograms();                             // fill from FlowEvent
00082     if (pFlowEvent) FillParticleHistograms();          // fill particle histograms
00083     if (Debug()) StMaker::PrintInfo();
00084   } else {
00085     gMessMgr->Info("##### FlowCumulantMaker: FlowEvent pointer null");
00086   }
00087   
00088   return kStOK;
00089 }
00090 
00091 //-----------------------------------------------------------------------
00092 
00093 Int_t StFlowCumulantMaker::Init() {
00094   // Book histograms
00095 
00096   float ptMaxPart = Flow::ptMaxPart;
00097   if (pFlowSelect->PtMaxPart()) {
00098     ptMaxPart = pFlowSelect->PtMaxPart();
00099   }
00100 
00101   nPtBinsPart = Flow::nPtBinsPart;
00102   if (pFlowSelect->PtBinsPart()) {
00103     nPtBinsPart = pFlowSelect->PtBinsPart();
00104   }
00105 
00106   xLabel = "Pseudorapidity";
00107   if (strlen(pFlowSelect->PidPart()) != 0) { xLabel = "Rapidity"; }  
00108   
00109   profScale=1000.;// in new root version, error does not show correctly 
00110   //if the number is too small. this is a trick to walkaround.
00111 
00112   r0 = 1.5; // this number should be small, but it could bring numerical 
00113                   // error if it is too small.
00114   r0Sq = r0 * r0;
00115   r0Mix= 3.;  //r0 for v1{3} calculation
00116 
00117   m_M = 1;  // if m_M = 2, what measured is v2{v1v1v2}, v4{v2v2v4}, v6{v3v3v6}... etc. 
00118   
00119   bool noDenomFileWarned = kFALSE;
00120   TString* histTitle;
00121 
00122   for (int k = 0; k < Flow::nSels; k++) {
00123 
00124     // for each selection
00125     histFull[k].mHistCumul = new TProfile*[Flow::nCumulDiffOrders];
00126     
00127     // mixed cumulant. (1st har mix with 2nd har. for v1 analysis).
00128     //we study 3-part v1 only, so we do not need sth. like nCumulDiffOrders as above.
00129     histTitle = new TString("Flow_CumulMix");
00130     histTitle->Append("_Sel");                      
00131     *histTitle +=k+1;
00132     histFull[k].mHistCumulMix =         
00133         new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5,
00134                      (float)(Flow::nHars) + 0.5, -1.*FLT_MAX, FLT_MAX, "");
00135     histFull[k].mHistCumulMix->SetXTitle("place for saving mixed cumulant");
00136     delete histTitle;    
00137 
00138 
00139     for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
00140       char theCumulOrder[2]; // if >10, need to use char*
00141       sprintf(theCumulOrder,"%d",(ord+1)*2);
00142       histTitle = new TString("Flow_Cumul_Order"); 
00143       *histTitle->Append(*theCumulOrder);           
00144       histTitle->Append("_Sel");                      
00145       *histTitle += k+1;          
00146       histFull[k].mHistCumul[ord] =  
00147         new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5,
00148                      (float)(Flow::nHars) + 0.5, -1.*FLT_MAX, FLT_MAX, "");
00149       histFull[k].mHistCumul[ord]->SetXTitle("harmonic");
00150       delete histTitle;
00151     }
00152     
00153     // for each harmonic
00154     for (int j = 0; j < Flow::nHars; j++) {
00155 
00156       histTitle = new TString("Flow_CumulMultSum_Sel");
00157       *histTitle += k+1;
00158       histTitle->Append("_Har");
00159       *histTitle += j+1;
00160       histFull[k].histFullHar[j].mHistMultSum =
00161         new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
00162       delete histTitle;
00163 
00164       histTitle = new TString("Flow_CumulMeanWgtSqrSum_Sel");
00165       *histTitle +=k+1;
00166       histTitle->Append("_Har");
00167       *histTitle +=j+1;
00168       histFull[k].histFullHar[j].mHistMeanWgtSqrSum =
00169         new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
00170       delete histTitle;
00171 
00172 
00173       histTitle = new TString("Flow_CumulNEvent_Sel");
00174       *histTitle += k+1;
00175       histTitle->Append("_Har");
00176       *histTitle += j+1;
00177       histFull[k].histFullHar[j].mHistNEvent =
00178         new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
00179       delete histTitle;
00180 
00181       // ****  for differential flow   ****
00182       
00183       // cumulant        d_p/n{k}      
00184       // p is the harmonic of the differential flow 
00185       // n is the harmonic of the integrated flow used for the differential flow
00186       // measurment. p = mn = either n or 2n.
00187       
00188       // if m_M=1, k=2    dn/n{2}  : cumulant from 2-part corr.
00189       // if m_M=1, k=4    dn/n{4}  : cumulant from 4-part corr.
00190       // if m_M=2, k=2    d2n/n{3} : cumulant from 3-part corr., mixed harmonic
00191       // if m_M=2, k=4    d2n/n{5} : cumulant from 5-part corr., mixed harmonic
00192       // where {2},{3} corresponds to theCumulIndex=1 below.
00193       // {4},{5} corresponds to theCumulIndex=2 below.
00194       
00195       histFull[k].histFullHar[j].mHistCumul2D  =  
00196         new TProfile2D*[Flow::nCumulDiffOrders];
00197       histFull[k].histFullHar[j].mHistCumulEta = 
00198         new TProfile*[Flow::nCumulDiffOrders];
00199       histFull[k].histFullHar[j].mHistCumulPt  = 
00200         new TProfile*[Flow::nCumulDiffOrders];
00201       
00202       // For each cumulant order
00203       for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
00204         char theCumulOrder[2]; // if >10, need to use char*
00205         sprintf(theCumulOrder,"%d",(ord+1)*2);
00206         
00207         histTitle = new TString("Flow_Cumul2D_Order");
00208         histTitle->Append(*theCumulOrder);
00209         histTitle->Append("_Sel");
00210         *histTitle += k+1;
00211         histTitle->Append("_Har");
00212         *histTitle += j+1;
00213         histFull[k].histFullHar[j].mHistCumul2D[ord] =  
00214           new TProfile2D(histTitle->Data(),histTitle->Data(), mNEtaBins,
00215                          mEtaMin, mEtaMax, nPtBinsPart, Flow::ptMin,
00216                          ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
00217         histFull[k].histFullHar[j].mHistCumul2D[ord]->SetXTitle((char*)xLabel.Data());
00218         histFull[k].histFullHar[j].mHistCumul2D[ord]->SetYTitle("Pt (GeV/c)");
00219         delete histTitle;
00220         
00221         histTitle = new TString("Flow_CumulEta_Order");
00222         histTitle->Append(*theCumulOrder);
00223         histTitle->Append("_Sel");
00224         *histTitle += k+1;
00225         histTitle->Append("_Har");
00226         *histTitle += j+1;
00227         histFull[k].histFullHar[j].mHistCumulEta[ord] =  
00228           new TProfile(histTitle->Data(),histTitle->Data(), mNEtaBins,
00229                        mEtaMin, mEtaMax, -1.*FLT_MAX, FLT_MAX, "");
00230         histFull[k].histFullHar[j].mHistCumulEta[ord]->SetXTitle((char*)xLabel.Data());
00231         delete histTitle;
00232         
00233         histTitle = new TString("Flow_CumulPt_Order");
00234         histTitle->Append(*theCumulOrder);
00235         histTitle->Append("_Sel");
00236         *histTitle += k+1;
00237         histTitle->Append("_Har");
00238         *histTitle += j+1;
00239         histFull[k].histFullHar[j].mHistCumulPt[ord] =  
00240           new TProfile(histTitle->Data(), histTitle->Data(), nPtBinsPart,
00241                        Flow::ptMin, ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
00242         histFull[k].histFullHar[j].mHistCumulPt[ord]->SetXTitle("Pt (GeV/c)");
00243         delete histTitle;
00244         
00245       }
00246       
00247 
00248       // for mixing 1st and 2nd harmonic. only coutHars=1 is meaningful.
00249       histTitle = new TString("Flow_CumulMix2D");
00250       histTitle->Append("_Sel");
00251       *histTitle +=k+1;
00252       histTitle->Append("_Har");
00253       *histTitle +=j+1;
00254       histFull[k].histFullHar[j].mHistCumulMix2D  =  
00255         new TProfile2D(histTitle->Data(),histTitle->Data(), Flow::nEtaBins,
00256                        Flow::etaMin, Flow::etaMax, nPtBinsPart, Flow::ptMin,
00257                        ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
00258       histFull[k].histFullHar[j].mHistCumulMix2D->SetXTitle((char*)xLabel.Data());
00259       histFull[k].histFullHar[j].mHistCumulMix2D->SetYTitle("Pt (GeV)");
00260       delete histTitle;
00261       
00262       
00263       histTitle = new TString("Flow_CumulMixEta");
00264       histTitle->Append("_Sel");
00265       *histTitle +=k+1;
00266       histTitle->Append("_Har");
00267       *histTitle +=j+1;
00268       histFull[k].histFullHar[j].mHistCumulMixEta =  
00269         new TProfile(histTitle->Data(),histTitle->Data(), Flow::nEtaBins,
00270                      Flow::etaMin, Flow::etaMax, -1.*FLT_MAX, FLT_MAX, "");
00271       histFull[k].histFullHar[j].mHistCumulMixEta->SetXTitle((char*)xLabel.Data());
00272       delete histTitle;
00273       
00274       
00275       histTitle = new TString("Flow_CumulMixPt");
00276       histTitle->Append("_Sel");
00277       *histTitle +=k+1;
00278       histTitle->Append("_Har");
00279       *histTitle +=j+1;
00280       histFull[k].histFullHar[j].mHistCumulMixPt =  
00281         new TProfile(histTitle->Data(), histTitle->Data(), nPtBinsPart,
00282                      Flow::ptMin, ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
00283       histFull[k].histFullHar[j].mHistCumulMixPt->SetXTitle("Pt (GeV)");
00284       delete histTitle;
00285       
00286       
00287       histFull[k].histFullHar[j].mCumulG0Denom = 
00288         new TProfile*[Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax];
00289 
00290       // For each cumulant order * qMax orders
00291       for (int pq  = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) { 
00292         
00293         TString* histTitleIntegDenom; //for read in
00294         
00295         int cumulIndex = (pq/Flow::nCumulDiff_qMax) + 1; // like 1,2,3.
00296         // for cumulant order 2,4,6 not begining with 0. This is "p" in Eq. (B1, PG5).
00297         int qIndex        = pq%Flow::nCumulDiff_qMax; // like 0,1,..._qMax-1  
00298         // begin with 0. This is "q" in Eq. (B3, PG5).
00299         
00300         char theCumulOrderChar[2];
00301         char qIndexOrderChar[2]; // if >10, need to use char*
00302         sprintf(theCumulOrderChar,"%d",cumulIndex*2); 
00303         sprintf(qIndexOrderChar,"%d",qIndex);
00304                 
00305         histTitle = new TString("Flow_CumulDenom_Order");
00306         histTitle->Append(*theCumulOrderChar);
00307         histTitle->Append("_GenFcnqIdx");
00308         histTitle->Append(*qIndexOrderChar);
00309         histTitle->Append("_Sel");
00310         *histTitle += k+1;
00311         histTitle->Append("_Har");
00312         *histTitle += j+1;
00313         histFull[k].histFullHar[j].mCumulG0Denom[pq] = 
00314           new TProfile(histTitle->Data(),histTitle->Data(), 1,0., 1., -1.*FLT_MAX, FLT_MAX, "");
00315         histFull[k].histFullHar[j].mCumulG0Denom[pq]->SetYTitle("<G>");
00316         histTitleIntegDenom = new TString(histTitle->Data());
00317         delete histTitle;
00318         
00319         // Open the file and get the histograms 
00320         // do not modify this section [ROOT bug] 
00321         TFile f("denominator.root","R");
00322         if (f.IsOpen()) {
00323           
00324           f.cd();
00325           TProfile* tempDenomProfile = 
00326             dynamic_cast<TProfile*>(f.Get(histTitleIntegDenom->Data()));
00327           if (!tempDenomProfile) {
00328             cout << "##### FlowCumulantAnalysis: can not find " <<
00329               histTitleIntegDenom->Data() << endl;
00330             return kFALSE;
00331           }   
00332           delete  histTitleIntegDenom;     
00333 
00334           histFull[k].histFullHar[j].mCumulG0DenomRead[pq]
00335             = tempDenomProfile->GetBinContent(1);
00336           
00337           f.Close();
00338           
00339         } else {
00340 
00341           if(!noDenomFileWarned) {
00342             gMessMgr->Info("##### FlowCumulantAnalysis:denominator.root is not present, assumming this run is just for producing denominator.root. That means cumulant flow result in flow.cumulant.root is nonsense for this run. ");
00343             noDenomFileWarned = kTRUE;
00344           }
00345           histFull[k].histFullHar[j].mCumulG0DenomRead[pq] = 1.; 
00346 
00347         }
00348         
00349         double theTempPhi = twopi*((double)qIndex) /  // Eq. (PG5)
00350           ((double)Flow::nCumulDiff_qMax); 
00351         double theRz = r0*::sqrt((double)cumulIndex);
00352         histFull[k].histFullHar[j].mDiffXz[pq] = theRz*cos(theTempPhi);
00353         histFull[k].histFullHar[j].mDiffYz[pq] = theRz*sin(theTempPhi);
00354       }
00355       
00356 
00357       histFull[k].histFullHar[j].mCumulG0MixDenom = 
00358         new   TProfile*[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
00359 
00360       histFull[k].histFullHar[j].mHistCumulIntegG0MixSum =
00361         new       TH1D*[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
00362 
00363 
00364       for (int pq  = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) { 
00365 
00366         
00367         TString* histTitleIntegDenomMix; //for read in
00368         
00369         int pIndex = (pq/Flow::nCumulMixHar_qMax);
00370         int qIndex        = pq%Flow::nCumulMixHar_qMax;
00371         
00372         char pIndexChar[2];
00373         char qIndexChar[2]; // if >10, need to use char*
00374         sprintf(pIndexChar,"%d",pIndex); 
00375         sprintf(qIndexChar,"%d",qIndex);
00376         
00377 
00378         histTitle = new TString("Flow_CumulMixDenom");
00379         histTitle->Append("_GenFunIdxp");
00380         histTitle->Append(*pIndexChar);
00381         histTitle->Append("_GenFunIdxq");
00382         histTitle->Append(*qIndexChar);
00383         histTitle->Append("_Sel");
00384         *histTitle +=k+1;
00385         histTitle->Append("_Har");
00386         *histTitle +=j+1;
00387         histFull[k].histFullHar[j].mCumulG0MixDenom[pq] =
00388           new TProfile(histTitle->Data(),histTitle->Data(), 1,0., 1., -1.*FLT_MAX, FLT_MAX, "");
00389         histFull[k].histFullHar[j].mCumulG0MixDenom[pq]->SetYTitle("<G>");
00390         histTitleIntegDenomMix = new TString(histTitle->Data());
00391         delete histTitle;
00392                                                                      
00393 
00394         histTitle = new TString("Flow_CumulIntegG0MixSum");
00395         histTitle->Append("_GenFunIdxp");
00396         histTitle->Append(*pIndexChar);
00397         histTitle->Append("_GenFunIdxq");
00398         histTitle->Append(*qIndexChar);
00399         histTitle->Append("_Sel");
00400         *histTitle +=k+1;
00401         histTitle->Append("_Har");
00402         *histTitle +=j+1;
00403         histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq] =
00404           new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
00405         delete histTitle;   
00406 
00407 
00408         histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] = 0.;
00409 
00410 
00411         TFile f("denominator.root","R");
00412         if (f.IsOpen()) {
00413           
00414           f.cd();
00415           TProfile* tempDenomMixProfile = 
00416             dynamic_cast<TProfile*>(f.Get(histTitleIntegDenomMix->Data()));
00417           if (!tempDenomMixProfile) {
00418             cout << "##### FlowCumulantAnalysis: can not find " <<
00419               histTitleIntegDenomMix->Data() << endl;
00420             return kFALSE;
00421           }   
00422           delete  histTitleIntegDenomMix;     
00423 
00424           histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq]
00425             = tempDenomMixProfile->GetBinContent(1);
00426           
00427           f.Close();
00428           
00429         } else {
00430 
00431           if(!noDenomFileWarned) {
00432             gMessMgr->Info("##### FlowCumulantAnalysis:denominator.root is not present, assumming this run is just for producing denominator.root. That means cumulant flow result with mixed harmonics in flow.cumulant.root is nonsense for this run. ");
00433             noDenomFileWarned = kTRUE;
00434           }
00435           histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq] = 1.; 
00436 
00437         }
00438                                                                                                                                     
00439     }
00440 
00441 
00442     for (int p = 0; p < Flow::nCumulMixHar_pMax; p++){ //v1{3} paper Eq (29)
00443       histFull[k].histFullHar[j].mMixX1z[p]=r0Mix*cos( pi*((double)p)/4. ); 
00444       histFull[k].histFullHar[j].mMixY1z[p]=r0Mix*sin( pi*((double)p)/4. ); 
00445     }
00446 
00447     for (int q = 0; q < Flow::nCumulMixHar_qMax; q++){
00448       histFull[k].histFullHar[j].mMixX2z[q]=r0Mix*cos( pi*((double)q)/2. ); 
00449       histFull[k].histFullHar[j].mMixY2z[q]=r0Mix*sin( pi*((double)q)/2. ); 
00450     }
00451 
00452 
00453     // ***  for integrated flow  ***
00454 
00455     histFull[k].histFullHar[j].mHistCumulIntegG0Sum =
00456       new TH1D*[Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax];
00457 
00458     for (int pq = 0; pq <  Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax; pq++) {
00459       int cumulIndex = (pq/Flow::nCumulInteg_qMax) + 1; // like 1,2,3. 
00460       //not begining with 0. That is "p" in Eq. (B1, PG5).
00461       int qIndex = pq%(Flow::nCumulInteg_qMax); // like 0,1,..qMax-1 
00462       //begining with 0. Just like Eq. (B3, PG5).
00463       
00464       char theCumulOrderChar[2];
00465       char qIndexOrderChar[2]; // if >10, need to use char*
00466       sprintf(theCumulOrderChar,"%d",cumulIndex*2); 
00467       sprintf(qIndexOrderChar,"%d",qIndex);
00468       
00469       double theTempPhi = twopi*((double)qIndex) /
00470           ((double)Flow::nCumulInteg_qMax); 
00471       double theRz = r0*::sqrt((double)cumulIndex);
00472       histFull[k].histFullHar[j].mIntegXz[pq] = theRz*cos(theTempPhi);
00473       histFull[k].histFullHar[j].mIntegYz[pq] = theRz*sin(theTempPhi);
00474       
00475       histFull[k].histFullHar[j].mCumulIntegG0[pq] = 0.;
00476 
00477       histTitle = new TString("Flow_CumulIntegG0Sum_Order");
00478       histTitle->Append(*theCumulOrderChar);
00479       histTitle->Append("_GenFunIdx");
00480       histTitle->Append(*qIndexOrderChar);
00481       histTitle->Append("_Sel");
00482       *histTitle += k+1;
00483       histTitle->Append("_Har");
00484       *histTitle += j+1;
00485       histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq] =
00486         new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
00487       delete histTitle;   
00488       
00489     }
00490     
00491     histFull[k].histFullHar[j].mMultSum       = 0.;
00492     histFull[k].histFullHar[j].mNEvent        = 0;
00493     histFull[k].histFullHar[j].mMeanWgtSqrSum = 0.;
00494     
00495     }
00496   }
00497   
00498   gMessMgr->SetLimit("##### FlowCumulantAnalysis", 2);
00499   gMessMgr->Info("##### FlowCumulantAnalysis: $Id: StFlowCumulantMaker.cxx,v 1.22 2006/02/22 19:12:29 posk Exp $");
00500 
00501   return StMaker::Init();
00502 }
00503 
00504 //-----------------------------------------------------------------------
00505 
00506 void StFlowCumulantMaker::FillFromFlowEvent() {
00507   // Get event quantities from StFlowEvent
00508   
00509   for (int k = 0; k < Flow::nSels; k++) {
00510     pFlowSelect->SetSelection(k);
00511     for (int j = 0; j < Flow::nHars; j++) {
00512       pFlowSelect->SetHarmonic(j);
00513       pFlowSelect->SetSubevent(-1);
00514       
00515       // full event quantities
00516       mMult[k][j]       = pFlowEvent->Mult(pFlowSelect);
00517       mSqrtOfSumWgtSqr[k][j] = ::sqrt(pFlowEvent->SumWeightSquare(pFlowSelect));
00518       
00519     }
00520   }
00521 }
00522 
00523 //-----------------------------------------------------------------------
00524 
00525 void StFlowCumulantMaker::FillEventHistograms() {
00526   // Fill histograms with event quantities
00527   
00528   for (int k = 0; k < Flow::nSels; k++) {
00529     pFlowSelect->SetSelection(k);
00530     for (int j = 0; j < Flow::nHars; j++) {
00531       pFlowSelect->SetHarmonic(j);
00532       
00533       histFull[k].histFullHar[j].mMultSum       += (float)mMult[k][j];
00534       histFull[k].histFullHar[j].mNEvent++;
00535       histFull[k].histFullHar[j].mMeanWgtSqrSum += 
00536          (mSqrtOfSumWgtSqr[k][j]*mSqrtOfSumWgtSqr[k][j])/(float)mMult[k][j];
00537 
00538       
00539       for (int pq = 0; pq < Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax; pq++) {   
00540           histFull[k].histFullHar[j].mCumulIntegG0[pq] += 
00541             pFlowEvent->G_New( pFlowSelect,  
00542                                histFull[k].histFullHar[j].mIntegXz[pq], 
00543                                histFull[k].histFullHar[j].mIntegYz[pq] );
00544       }
00545 
00546       for (int pq  = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) { 
00547         int pIndex = pq/Flow::nCumulMixHar_qMax;
00548         int qIndex = pq%Flow::nCumulMixHar_qMax;
00549          if (j ==0)       //only for v1
00550          histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] +=
00551             pFlowEvent->G_Mix( pFlowSelect,        
00552                                histFull[k].histFullHar[j].mMixX1z[pIndex],
00553                                histFull[k].histFullHar[j].mMixY1z[pIndex],
00554                                histFull[k].histFullHar[j].mMixX2z[qIndex],
00555                                histFull[k].histFullHar[j].mMixY2z[qIndex] );
00556      }
00557 
00558     }
00559   }
00560   
00561 }
00562 
00563 //-----------------------------------------------------------------------
00564 
00565 void StFlowCumulantMaker::FillParticleHistograms() {
00566   // Fill histograms from the particles
00567   
00568   // Initialize Iterator
00569   StFlowTrackCollection* pFlowTracks = pFlowEvent->TrackCollection();
00570   StFlowTrackIterator itr;
00571   
00572   // In the view of good coding, the following block should be placed inside
00573   // the track loop. It is put here to reduce run time.
00574   // If somebody changes the pFlowSelect to select tracks for v-part, 
00575   // he/she has to change pFlowSelect in this block for consistency. 
00576   // Selections have to be consistent between here and the loop.
00577 
00578   double* evtG[Flow::nSels][Flow::nHars];
00579   double* theCrossterm[Flow::nSels][Flow::nHars];
00580   
00581   double* evtGMix[Flow::nSels][Flow::nHars];
00582   double* theCrosstermMix[Flow::nSels][Flow::nHars];
00583 
00584   for (int k = 0; k < Flow::nSels; k++) {
00585     for (int j = 0; j < Flow::nHars; j++) { 
00586        evtG[k][j] = 
00587         new double[Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax];
00588        theCrossterm[k][j]    = 
00589         new double[Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax];
00590        evtGMix[k][j] = 
00591         new double[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
00592        theCrosstermMix[k][j]    = 
00593         new double[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
00594     }
00595   }
00596 
00597   for (int k = 0; k < Flow::nSels; k++) { 
00598     for (int j = 0; j < Flow::nHars; j++) {
00599       pFlowSelect->SetSelection(k);
00600       pFlowSelect->SetHarmonic(j);
00601       
00602       
00603       for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
00604         
00605         evtG[k][j][pq] = 
00606           (pFlowEvent->G_New( pFlowSelect,
00607                               histFull[k].histFullHar[j].mDiffXz[pq],
00608                               histFull[k].histFullHar[j].mDiffYz[pq] )) ;
00609         theCrossterm[k][j][pq] =evtG[k][j][pq];
00610         histFull[k].histFullHar[j].mCumulG0Denom[pq]->Fill(0.5,evtG[k][j][pq]);
00611       }
00612 
00613       for (int pq  = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) { 
00614         int pIndex = pq/Flow::nCumulMixHar_qMax;
00615         int qIndex = pq%Flow::nCumulMixHar_qMax;
00616         if (j ==0) //only for v1          
00617         evtGMix[k][j][pq] =
00618             pFlowEvent->G_Mix( pFlowSelect,        
00619                                histFull[k].histFullHar[j].mMixX1z[pIndex],
00620                                histFull[k].histFullHar[j].mMixY1z[pIndex],
00621                                histFull[k].histFullHar[j].mMixX2z[qIndex],
00622                                histFull[k].histFullHar[j].mMixY2z[qIndex] );
00623         theCrosstermMix[k][j][pq] =evtGMix[k][j][pq];
00624         histFull[k].histFullHar[j].mCumulG0MixDenom[pq]->Fill(0.5,evtGMix[k][j][pq]);
00625       }  
00626  
00627     }
00628   }
00629 
00630   for (itr = pFlowTracks->begin(); itr != pFlowTracks->end(); itr++) {
00631     StFlowTrack* pFlowTrack = *itr;
00632     
00633     float phi = pFlowTrack->Phi();
00634     if (phi < 0.) phi += twopi;
00635     float eta = pFlowTrack->Eta();
00636     float pt  = pFlowTrack->Pt();
00637     
00638     for (int k = 0; k < Flow::nSels; k++) {
00639       pFlowSelect->SetSelection(k);
00640       double cumuTemp[Flow::nCumulDiffOrders];
00641       double cumuTempFlip[Flow::nCumulDiffOrders];
00642       double cumuTempMix;
00643       double cumuTempMixFlip;
00644       double order = 0.;
00645       double phiWgt = 1.;
00646       double phiWgtRaw = 1.; //no need raw phi weight for cumulants.
00647 
00648       for (int j = 0; j < Flow::nHars; j++) {
00649         bool oddHar = (j+1) % 2;
00650         pFlowSelect->SetHarmonic(j);
00651         order  = (double)(j+1);
00652         
00653         if (pFlowSelect->Select(pFlowTrack)) {
00654                   // Get phiWgt
00655           phiWgt = pFlowEvent->Weight(k, j, pFlowTrack);
00656         }
00657         
00658         // Caculate v for all particles selected for correlation analysis
00659         if (pFlowSelect->SelectPart(pFlowTrack)) {
00660           
00661           float yOrEta = 
00662             (strlen(pFlowSelect->PidPart()) != 0) ?  pFlowTrack->Y() : eta;
00663           
00664           double Dp[Flow::nCumulDiffOrders]; // the Dp in (B7, PG11)
00665           for (int pq = 0; pq < Flow::nCumulDiffOrders; pq++) {
00666             Dp[pq] = 0.; }
00667           
00668           for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
00669             int theCumulOrder = (pq/Flow::nCumulDiff_qMax) + 1; // like 1,2. 
00670             // not begin with 0. which is "p" in (B6, PG5)
00671             int qIndex        = pq%(Flow::nCumulDiff_qMax); // like 0,1,..qMax-1 
00672             // begin with 0. which is "q" in (B6, PG5)
00673             
00674             //Calculate PG11
00675             double theCoeff = ::pow(r0*::sqrt(double(theCumulOrder)), double(m_M)) /
00676               float(Flow::nCumulDiff_qMax); // first term in (B7)
00677             double theCosTerm = cos(twopi*float(qIndex)*float(m_M) /
00678                                     float(Flow::nCumulDiff_qMax)); // cos() in (B7)
00679             double theSinTerm = sin(twopi*float(qIndex)*float(m_M) /
00680                                     float(Flow::nCumulDiff_qMax)); // sin() in (B7)
00681                     
00682             if ( (pFlowSelect->SelectPart(pFlowTrack)) && 
00683                  (pFlowSelect->Select(pFlowTrack)) ) { // remove autocorrelation
00684 
00685                 theCrossterm[k][j][pq] = evtG[k][j][pq] / 
00686                   (1. + (phiWgt/mMult[k][j]) *
00687                    (2.*histFull[k].histFullHar[j].mDiffXz[pq]*cos(phi*order) +
00688                     2.*histFull[k].histFullHar[j].mDiffYz[pq]*sin(phi * order) ) ); 
00689                 // the argument in the last paragraph in page 9. 
00690               
00691             }
00692 
00693             
00694             double theXpq = (theCrossterm[k][j][pq]*cos(float(m_M) * order * phi)) / // (B6)
00695               histFull[k].histFullHar[j].mCumulG0DenomRead[pq]; 
00696             double theYpq = (theCrossterm[k][j][pq]*sin(float(m_M) * order * phi)) / // (B6)
00697               histFull[k].histFullHar[j].mCumulG0DenomRead[pq];
00698             
00699             Dp[theCumulOrder-1] += 
00700               theCoeff*(theCosTerm*theXpq + theSinTerm*theYpq); // (B7, PG11)
00701             
00702           }
00703           
00704 
00706 
00707           if (j==0){ //Har 1 only
00708 
00709           double DpxMix[Flow::nCumulMixHar_pMax]; 
00710           double DpyMix[Flow::nCumulMixHar_pMax]; 
00711 
00712           double DpqxMix[Flow::nCumulMixHar_pMax][Flow::nCumulMixHar_qMax]; 
00713           double DpqyMix[Flow::nCumulMixHar_pMax][Flow::nCumulMixHar_qMax]; 
00714 
00715 
00716           for (int p = 0; p < Flow::nCumulMixHar_pMax; p++) {
00717             DpxMix[p] = 0.;  DpyMix[p] = 0.;}
00718           
00719           for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax; pq++) {
00720 
00721             int pIndex = pq/Flow::nCumulMixHar_qMax;
00722             int qIndex = pq%Flow::nCumulMixHar_qMax;
00723 
00724             if ( (pFlowSelect->SelectPart(pFlowTrack)) && 
00725                  (pFlowSelect->Select(pFlowTrack)) ) { // rm autocorrelation
00726 
00727               double etaWgt = (oddHar) ? ( (eta>0) ? (pFlowEvent->EtaAbsWgtValue(eta)) :  (-1.*(pFlowEvent->EtaAbsWgtValue(eta))) ) : 1.;
00728 
00729               double ptWgt  = pFlowEvent->PtAbsWgtValue(pt);
00730               
00731               double detectorV1Wgt = 1.;
00732 
00733               if (pFlowTrack->TopologyMap().hasHitInDetector(kTpcId) ||  
00734                   (pFlowTrack->TopologyMap().data(0) == 0 && 
00735                    pFlowTrack->TopologyMap().data(1) == 0)) {
00736                 // hits in Tpc or TopologyMap not available
00737                 detectorV1Wgt =pFlowEvent->V1TPCDetctWgtG_Mix(k);
00738               } else if (pFlowTrack->TopologyMap().trackFtpcEast() ) {
00739                 detectorV1Wgt =pFlowEvent->V1FtpcEastDetctWgtG_Mix(k);
00740               } else if (pFlowTrack->TopologyMap().trackFtpcWest() ) { 
00741                 detectorV1Wgt =pFlowEvent->V1FtpcWestDetctWgtG_Mix(k);
00742               }
00743 
00744 
00745               double detectorV2Wgt = 1.;
00746 
00747               if (pFlowTrack->TopologyMap().hasHitInDetector(kTpcId) ||  
00748                   (pFlowTrack->TopologyMap().data(0) == 0 && 
00749                    pFlowTrack->TopologyMap().data(1) == 0)) {
00750                 // hits in Tpc or TopologyMap not available
00751                 detectorV2Wgt =pFlowEvent->V2TPCDetctWgtG_Mix(k);
00752               } else if (pFlowTrack->TopologyMap().trackFtpcEast() ) {
00753                 detectorV2Wgt =pFlowEvent->V2FtpcEastDetctWgtG_Mix(k);
00754               } else if (pFlowTrack->TopologyMap().trackFtpcWest() ) { 
00755                 detectorV2Wgt =pFlowEvent->V2FtpcWestDetctWgtG_Mix(k);
00756               }
00757               
00758 
00759               theCrosstermMix[k][j][pq] = evtGMix[k][j][pq] / 
00760                 (1. + (phiWgtRaw*etaWgt*detectorV1Wgt/mMult[k][j]) * 
00761                  (2.* histFull[k].histFullHar[j].mMixX1z[pIndex] * cos(phi * order) + 
00762                   2.*  histFull[k].histFullHar[j].mMixY1z[pIndex] * sin(phi * order) ) 
00763                  + (phiWgtRaw*ptWgt*detectorV2Wgt/mMult[k][j]) * 
00764                  (2.* histFull[k].histFullHar[j].mMixX2z[qIndex] * cos(phi * order*2.) + 
00765                   2.* histFull[k].histFullHar[j].mMixY2z[qIndex] * sin(phi * order*2.) ) );
00766               
00767             }
00768             
00769             // for writting out <G(p,q)>, the denominator 
00770             
00771             DpqxMix[pIndex][qIndex] = theCrosstermMix[k][j][pq] * cos(order * phi)/
00772               histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq];
00773             DpqyMix[pIndex][qIndex] = theCrosstermMix[k][j][pq] * sin(order * phi)/
00774               histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq];
00775             
00776           }
00777 
00778           
00779           for (int p=0; p< Flow::nCumulMixHar_pMax; p++){
00780             DpxMix[p]=(1./(4.*r0Mix))*( DpqxMix[p][0]-DpqxMix[p][2]+DpqyMix[p][1]-DpqyMix[p][3]);
00781             DpyMix[p]=(1./(4.*r0Mix))*( DpqyMix[p][0]-DpqyMix[p][2]+DpqxMix[p][3]-DpqxMix[p][1]);
00782           } 
00783           
00784 
00785           cumuTempMix = (1./(4.*r0Mix))*(DpxMix[0]-DpyMix[2]-DpxMix[4]+DpyMix[6]);
00786           
00787           cumuTempMixFlip = cumuTempMix;
00788           
00789           if (eta < 0 && oddHar) {
00790             cumuTempMixFlip *= -1.;
00791           }
00792           
00793 
00794           histFull[k].histFullHar[j].mHistCumulMix2D->
00795             Fill(yOrEta, pt, cumuTempMix*profScale); 
00796           histFull[k].histFullHar[j].mHistCumulMixEta->
00797             Fill(yOrEta, cumuTempMix*profScale); 
00798           histFull[k].histFullHar[j].mHistCumulMixPt->
00799             Fill(pt, cumuTempMixFlip*profScale); 
00800           
00801           histFull[k].mHistCumulMix->Fill(order,cumuTempMixFlip*profScale); //only order ==1 get filled anyway.
00802           
00803           }
00804 
00806 
00807 
00808           if (m_M==1) {
00809             cumuTemp[0] = ((2.*Dp[1-1])-(0.5*Dp[2-1]))/r0Sq; // (B9, PG12)
00810             cumuTemp[1] = ((-2.*Dp[1-1])+Dp[2-1])/(r0Sq*r0Sq);
00811           } else if (m_M==2) {
00812             cumuTemp[0] = ((4.*Dp[1-1])-(0.5*Dp[2-1]))/(r0Sq*r0Sq); // (B10)
00813             cumuTemp[1] = ((-6.*Dp[1-1])+(1.5*Dp[2-1]))/(r0Sq*r0Sq*r0Sq);
00814           }
00815           
00816           cumuTempFlip[0] = cumuTemp[0];
00817           cumuTempFlip[1] = cumuTemp[1];
00818           if (eta < 0 && oddHar) {
00819             cumuTempFlip[0] *= -1.;
00820             cumuTempFlip[1] *= -1.;
00821           }
00822           
00823           for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
00824             histFull[k].histFullHar[j].mHistCumul2D[ord]->
00825               Fill(yOrEta, pt,  ((ord>0) ? cumuTemp[ord]*profScale : cumuTemp[ord])); 
00826             histFull[k].histFullHar[j].mHistCumulEta[ord]->
00827               Fill(yOrEta, ((ord>0) ? cumuTemp[ord]*profScale : cumuTemp[ord])); 
00828             histFull[k].histFullHar[j].mHistCumulPt[ord]->
00829               Fill(pt, ((ord>0) ? cumuTempFlip[ord]*profScale : cumuTempFlip[ord])); 
00830           }       
00831           
00832           for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++)
00833             histFull[k].mHistCumul[ord]->Fill(order, ((ord>0) ? cumuTempFlip[ord]*profScale : cumuTempFlip[ord]) );
00834         }
00835       }
00836     }  
00837   }
00838   
00839 }
00840 
00841 //-----------------------------------------------------------------------
00842 
00843 Int_t StFlowCumulantMaker::Finish() {
00844   
00845   TString* histTitle;
00846   
00847   TOrdCollection* XpqYpqDenomNames = new TOrdCollection(Flow::nSels*Flow::nHars); 
00848   TOrdCollection* savedHistNames = new TOrdCollection(Flow::nSels*Flow::nHars*Flow::nCumulDiffOrders); 
00849   
00850   cout << endl << "##### Cumulant Maker:" << endl;
00851   for (int k = 0; k < Flow::nSels; k++) {
00852     
00853     cout << "##### selection "<< k+1 <<"  #### "<<endl;
00854     
00855     // integrated flow from cumulant
00856     histFull[k].mHist_v = new TH1D*[Flow::nCumulDiffOrders];
00857     
00858     for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
00859       char theCumulOrder[2]; // if >10, need to use char*
00860       sprintf(theCumulOrder,"%d",(ord+1)*2);
00861 
00862       histTitle = new TString("Flow_Cumul_v_Order");
00863       histTitle->Append(*theCumulOrder);
00864       histTitle->Append("_Sel");
00865       *histTitle += k+1;
00866       histFull[k].mHist_v[ord] = 
00867         new  TH1D(*(histFull[k].mHistCumul[ord]->ProjectionX(histTitle->Data(),"e")));
00868       histFull[k].mHist_v[ord]->SetTitle(histTitle->Data());
00869       histFull[k].mHist_v[ord]->SetXTitle("harmonic");
00870       histFull[k].mHist_v[ord]->SetYTitle("v (%)");
00871       delete histTitle;
00872       if (ord>0)  histFull[k].mHist_v[ord]->Scale(1./profScale);
00873       savedHistNames->AddLast(histFull[k].mHist_v[ord]);
00874     }
00875 
00877     histTitle = new TString("Flow_CumulMix_v_Sel");
00878     *histTitle +=k+1;   
00879     histFull[k].mHistMix_v = 
00880       new  TH1D(*(histFull[k].mHistCumulMix->ProjectionX(histTitle->Data(),"e")));
00881     histFull[k].mHistMix_v->SetTitle(histTitle->Data());
00882     histFull[k].mHistMix_v->SetXTitle("place for v1 from mixed harmonic");
00883     histFull[k].mHistMix_v->SetYTitle("v (%)");
00884     delete histTitle;
00885     
00886     histFull[k].mHistMix_v->Scale(1./profScale);
00887     savedHistNames->AddLast(histFull[k].mHistMix_v);
00888     
00890     
00891     double  meanIntegV[Flow::nHars];     // V**1
00892     double  meanIntegV2[Flow::nHars];    // V**2
00893     double  meanIntegV3[Flow::nHars];    // V**3
00894     double  meanIntegV4[Flow::nHars];    // V**4
00895     double  cumulInteg1[Flow::nHars];    // outside of harmonic loop
00896     double  cumulInteg2[Flow::nHars];
00897     double  cumulInteg3[Flow::nHars];
00898     
00899     double cumulIntegMix[Flow::nHars];
00900     double meanIntegVMix[Flow::nHars];
00901 
00902     for (int j = 0; j < Flow::nHars; j++) {
00903       meanIntegV[j]  = 0.;
00904       meanIntegV2[j] = 0.;
00905       meanIntegV3[j] = 0.;
00906       meanIntegV4[j] = 0.;
00907       cumulInteg1[j] = 0.;
00908       cumulInteg2[j] = 0.;
00909       cumulInteg3[j] = 0.;
00910       cumulIntegMix[j] = 0.;
00911       meanIntegVMix[j] = 0.;
00912     }
00913     
00914     for (int j = 0; j < Flow::nHars; j++) {
00915       
00916       double mAvMult = // average multiplicity
00917         float(histFull[k].histFullHar[j].mMultSum)/
00918         (float(histFull[k].histFullHar[j].mNEvent));
00919       
00920       histFull[k].histFullHar[j].mHistMultSum->
00921         SetBinContent(1,double(histFull[k].histFullHar[j].mMultSum));
00922       histFull[k].histFullHar[j].mHistNEvent->
00923         SetBinContent(1,double(histFull[k].histFullHar[j].mNEvent));
00924       histFull[k].histFullHar[j].mHistMeanWgtSqrSum->
00925         SetBinContent(1,double(histFull[k].histFullHar[j].mMeanWgtSqrSum));
00926       
00927       double CpInteg[Flow::nCumulIntegOrders]; // Cp in (B4, PG6)
00928       
00929       for (int pq = 0; pq < Flow::nCumulIntegOrders; pq ++) CpInteg[pq] = 0.;
00930       for (int pq = 0; pq < Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax; pq++) {   
00931         int theCumulOrder = (pq/Flow::nCumulInteg_qMax) + 1; // like 1,2,3. 
00932         // not begining with 0. That is "p" in (B3, PG5)
00933         //int qIndex = pq%(Flow::nCumulInteg_qMax); // like 0,1,..qMax-1 
00934         // begining with 0. "q" in (B3, PG5)
00935 
00936         histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq]-> // don't move it.
00937           SetBinContent(1,histFull[k].histFullHar[j].mCumulIntegG0[pq]);
00938 
00939         //mCumulIntegG0 is the sum of IntegG0 until the code right below,
00940         //where mCumulIntegG0 becomes <Gn(z)> :
00941         histFull[k].histFullHar[j].mCumulIntegG0[pq] /= 
00942           float(histFull[k].histFullHar[j].mNEvent);        // <Gn(z)> (PG 4)
00943 
00944         CpInteg[theCumulOrder-1] +=
00945           (mAvMult*(::pow(histFull[k].histFullHar[j].mCumulIntegG0[pq], 1./mAvMult) -1.) /
00946            float(Flow::nCumulInteg_qMax)); // (B3, PG6) 
00947         
00948       }
00949       
00950       // add Xpq Ypq denominator to write out file list
00951       for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
00952         XpqYpqDenomNames->AddLast(histFull[k].histFullHar[j].mCumulG0Denom[pq]);
00953       }
00954 
00955       cumulInteg1[j] = // (B5, PG7)
00956         (3.*CpInteg[1-1] - (3./2.)*CpInteg[2-1] + (1./3.)*CpInteg[3-1]) / r0Sq;
00957       
00958       cumulInteg2[j] = ((-10.*CpInteg[1-1]) + (8.*CpInteg[2-1]) - 
00959                          (2.*CpInteg[3-1])) / (r0Sq*r0Sq); 
00960       
00961       cumulInteg3[j] = ( (18.*CpInteg[1-1]) - (18.*CpInteg[2-1]) + (6.*CpInteg[3-1]))
00962                         / (r0Sq*r0Sq*r0Sq);
00963       
00964       // now histograms for flow results:
00965       histFull[k].histFullHar[j].mHist_v2D  = new TH2D*[Flow::nCumulDiffOrders];
00966       histFull[k].histFullHar[j].mHist_vEta = new TH1D*[Flow::nCumulDiffOrders];
00967       histFull[k].histFullHar[j].mHist_vPt  = new TH1D*[Flow::nCumulDiffOrders];
00968       
00969       for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
00970         char theCumulOrder[2]; // if >10, need to use char*
00971         sprintf(theCumulOrder,"%d",(ord+1)*2);
00972         
00973         histTitle = new TString("Flow_Cumul_v2D_Order");
00974         histTitle->Append(*theCumulOrder);
00975         histTitle->Append("_Sel");
00976         *histTitle += k+1;
00977         histTitle->Append("_Har");
00978         *histTitle += j+1;
00979         histFull[k].histFullHar[j].mHist_v2D[ord] = 
00980           new TH2D(*(histFull[k].histFullHar[j].mHistCumul2D[ord]->
00981                      ProjectionXY(histTitle->Data(),"e")));
00982         histFull[k].histFullHar[j].mHist_v2D[ord]->SetTitle(histTitle->Data());
00983         histFull[k].histFullHar[j].mHist_v2D[ord]->SetXTitle((char*)xLabel.Data());
00984         histFull[k].histFullHar[j].mHist_v2D[ord]->SetYTitle("Pt (GeV/c)");
00985         histFull[k].histFullHar[j].mHist_v2D[ord]->SetZTitle("v (%)");
00986         delete histTitle;
00987         
00988         if (ord>0) histFull[k].histFullHar[j].mHist_v2D[ord]->Scale(1./profScale);
00989 
00990         histTitle = new TString("Flow_Cumul_vEta_Order");
00991         histTitle->Append(*theCumulOrder);
00992         histTitle->Append("_Sel");
00993         *histTitle += k+1;
00994         histTitle->Append("_Har");
00995         *histTitle += j+1;
00996         histFull[k].histFullHar[j].mHist_vEta[ord] = 
00997           new  TH1D(*(histFull[k].histFullHar[j].mHistCumulEta[ord]->
00998                       ProjectionX(histTitle->Data(),"e")));
00999         histFull[k].histFullHar[j].mHist_vEta[ord]->SetTitle(histTitle->Data());
01000         histFull[k].histFullHar[j].mHist_vEta[ord]->SetXTitle((char*)xLabel.Data());
01001         histFull[k].histFullHar[j].mHist_vEta[ord]->SetYTitle("v (%)");
01002         delete histTitle;
01003         
01004         if (ord>0) histFull[k].histFullHar[j].mHist_vEta[ord]->Scale(1./profScale);
01005 
01006         histTitle = new TString("Flow_Cumul_vPt_Order");
01007         histTitle->Append(*theCumulOrder);
01008         histTitle->Append("_Sel");
01009         *histTitle += k+1;
01010         histTitle->Append("_Har");
01011         *histTitle += j+1;
01012         histFull[k].histFullHar[j].mHist_vPt[ord] = 
01013           new  TH1D(*(histFull[k].histFullHar[j].mHistCumulPt[ord]->
01014                       ProjectionX(histTitle->Data(),"e")));
01015         histFull[k].histFullHar[j].mHist_vPt[ord]->SetTitle(histTitle->Data());
01016         histFull[k].histFullHar[j].mHist_vPt[ord]->SetXTitle("Pt (GeV/c)");
01017         histFull[k].histFullHar[j].mHist_vPt[ord]->SetYTitle("v (%)");
01018         delete histTitle;
01019 
01020         if (ord>0) histFull[k].histFullHar[j].mHist_vPt[ord]->Scale(1./profScale);
01021 
01022       }
01023       
01024       // new method, Eq. (PG8)
01025       meanIntegV[j]  = ::sqrt(cumulInteg1[j]);           // <v>    2-part, m=1
01026       meanIntegV2[j] = cumulInteg1[j];                 // <v**2> 2-part, m=2
01027       meanIntegV3[j] = ::pow(-1.*cumulInteg2[j], 3./4.); // <v**3> 4-part, m=1
01028       meanIntegV4[j] = -1.*cumulInteg2[j];             // <v**4> 4-part, m=2
01029       
01030       if (meanIntegV2[j]<0.) cout<<" Sel"<<k+1<<", <V**2> less than zero ! v"
01031                                  <<j+1<<" from 2 particle correlation failed."<<endl;
01032       if (meanIntegV4[j]<0.) cout<<" Sel"<<k+1<<", <V**4> less than zero ! v"
01033                                  <<j+1<<" from 4 particle correlation failed."<<endl;
01034   
01035       
01036       if (m_M==1) { // Eq. (PG14)
01037         histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
01038         histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-1./(meanIntegV3[j]*perCent));// (34b)
01039         histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
01040         histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
01041         histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
01042         histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
01043       } else if (m_M==2) {
01044         histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
01045         histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
01046         histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV2[j]*perCent));   // (35a)
01047         histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-0.5/(meanIntegV4[j]*perCent) ); // (35b)
01048         histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
01049         histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
01050       }
01051       
01052       for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
01053         savedHistNames->AddLast(histFull[k].histFullHar[j].mHist_v2D[ord]);
01054         savedHistNames->AddLast(histFull[k].histFullHar[j].mHist_vEta[ord]);
01055         savedHistNames->AddLast(histFull[k].histFullHar[j].mHist_vPt[ord]);
01056       }
01057     }
01058 
01059 
01061 
01062     for (int j = 0; j < Flow::nHars; j++) {
01063       if (j != 0) continue; //only j==0 makes sense for mixed har.
01064           
01065       double mAvMult = // average multiplicity
01066         float(histFull[k].histFullHar[j].mMultSum)/
01067         (float(histFull[k].histFullHar[j].mNEvent));
01068 
01069       double mAveMeanWgtSqr = // <w**2>
01070         float(histFull[k].histFullHar[j].mMeanWgtSqrSum)/
01071         (float(histFull[k].histFullHar[j].mNEvent));
01072 
01073       double  CpqMix[Flow::nCumulMixHar_pMax][Flow::nCumulMixHar_qMax];//(30)
01074       for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax; pq++) {
01075         int pIndex = pq/Flow::nCumulMixHar_qMax;
01076         int qIndex = pq%Flow::nCumulMixHar_qMax;
01077         
01078         histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq]->
01079           SetBinContent(1,histFull[k].histFullHar[j].mCumulIntegG0Mix[pq]);
01080         histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] /=  //<G(z1,z2)>
01081           float(histFull[k].histFullHar[j].mNEvent);    
01082         
01083         CpqMix[pIndex][qIndex]=mAvMult*(pow(histFull[k].histFullHar[j].mCumulIntegG0Mix[pq], 1./mAvMult) -1.); //Mix (21)
01084         
01085       }
01086       
01087       double CpxMix[Flow::nCumulMixHar_pMax];
01088       double CpyMix[Flow::nCumulMixHar_pMax];
01089       
01090       for (int p =0; p<Flow::nCumulMixHar_pMax; p++){//Mix (31)
01091         CpxMix[p] = (1./(4.*r0Mix))*(CpqMix[p][0] - CpqMix[p][2]);
01092         CpyMix[p] = (1./(4.*r0Mix))*(CpqMix[p][3] - CpqMix[p][1]);
01093       }
01094 
01095       cumulIntegMix[j] = (1./(4.*r0Mix*r0Mix))*( // Mix (32)
01096                                                 CpxMix[0]-CpyMix[1]-CpxMix[2]+CpyMix[3]
01097                                                 +CpxMix[4]-CpyMix[5]-CpxMix[6]+CpyMix[7]);
01098 
01099  
01100       double  tempMeanV = pow(-1.*cumulInteg2[1]*mAveMeanWgtSqr*mAveMeanWgtSqr,1./4.); // <wgt*v2>
01101       double  tempMeanVMixSq = cumulIntegMix[j]/tempMeanV; //(24)
01102       
01103       if (tempMeanVMixSq>0.)
01104         meanIntegVMix[j] = sqrt(tempMeanVMixSq);//(24)
01105       else cout<<"### <wgt*v1>**2 = "<<tempMeanVMixSq<<" < 0. failed "<<endl;
01106 
01107 
01108       histTitle = new TString("Flow_CumulMix_v2D_Sel");
01109       *histTitle +=k+1;
01110       histTitle->Append("_Har");
01111       *histTitle +=j+1;
01112       histFull[k].histFullHar[j].mHistMix_v2D = 
01113         new TH2D(*(histFull[k].histFullHar[j].mHistCumulMix2D->
01114                    ProjectionXY(histTitle->Data(),"e")));
01115       histFull[k].histFullHar[j].mHistMix_v2D->SetTitle(histTitle->Data());
01116       histFull[k].histFullHar[j].mHistMix_v2D->SetXTitle((char*)xLabel.Data());
01117       histFull[k].histFullHar[j].mHistMix_v2D->SetYTitle("Pt (GeV)");
01118       histFull[k].histFullHar[j].mHistMix_v2D->SetZTitle("v (%)");
01119       delete histTitle;
01120       
01121       histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./profScale);
01122       
01123       histTitle = new TString("Flow_CumulMix_vEta_Sel");
01124       *histTitle +=k+1;
01125       histTitle->Append("_Har");
01126       *histTitle +=j+1;
01127       histFull[k].histFullHar[j].mHistMix_vEta = 
01128         new  TH1D(*(histFull[k].histFullHar[j].mHistCumulMixEta->
01129                     ProjectionX(histTitle->Data(),"e")));
01130       histFull[k].histFullHar[j].mHistMix_vEta->SetTitle(histTitle->Data());
01131       histFull[k].histFullHar[j].mHistMix_vEta->SetXTitle((char*)xLabel.Data());
01132       histFull[k].histFullHar[j].mHistMix_vEta->SetYTitle("v (%)");
01133       delete histTitle;
01134       
01135       histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./profScale);
01136       
01137 
01138       histTitle = new TString("Flow_CumulMix_vPt_Sel");
01139       *histTitle +=k+1;
01140       histTitle->Append("_Har");
01141       *histTitle +=j+1;
01142       histFull[k].histFullHar[j].mHistMix_vPt = 
01143         new  TH1D(*(histFull[k].histFullHar[j].mHistCumulMixPt->
01144                     ProjectionX(histTitle->Data(),"e")));
01145       histFull[k].histFullHar[j].mHistMix_vPt->SetTitle(histTitle->Data());
01146       histFull[k].histFullHar[j].mHistMix_vPt->SetXTitle("Pt (GeV)");
01147       histFull[k].histFullHar[j].mHistMix_vPt->SetYTitle("v (%)");
01148       delete histTitle;
01149       
01150       histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./profScale);
01151       
01152       histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent)); 
01153       histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
01154       histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent)); 
01155       histFull[k].mHistMix_v->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent)); 
01156       
01157       savedHistNames->AddLast(histFull[k].histFullHar[j].mHistMix_v2D);
01158       savedHistNames->AddLast(histFull[k].histFullHar[j].mHistMix_vEta);
01159       savedHistNames->AddLast(histFull[k].histFullHar[j].mHistMix_vPt);
01160       savedHistNames->AddLast(histFull[k].mHistMix_v);
01161       
01162     } 
01164     
01165     
01166     if (m_M==1) {
01167       
01168       TH1D* histOfMeanIntegV = new TH1D(*(histFull[k].mHist_v[0]));
01169       histOfMeanIntegV->Reset();
01170       
01171       TH1D* histOfMeanIntegV3 = new TH1D(*(histFull[k].mHist_v[1]));
01172       histOfMeanIntegV3->Reset();
01173       
01174       for (int j = 1; j < Flow::nHars+1; j++) {
01175         histOfMeanIntegV->SetBinContent(j, 1./(meanIntegV[j-1]*perCent));
01176         histOfMeanIntegV->SetBinError(j,0.);
01177         histOfMeanIntegV3->SetBinContent(j, -1./(meanIntegV3[j-1]*perCent));
01178         histOfMeanIntegV3->SetBinError(j,0.);
01179       }
01180       histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV);
01181       histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV3);
01182       
01183       //force to be zero if the value is "nan"
01184       for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++){
01185         for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
01186           if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
01187                  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ) {
01188             histFull[k].mHist_v[ord]->SetBinContent(j,0.);
01189             histFull[k].mHist_v[ord]->SetBinError(j,0.);
01190           }
01191         }
01192       }
01193 
01194       for (int j = 1; j < Flow::nHars+1; j++) {
01195         cout << "##### 2-part v" << j << " = (" 
01196              << histFull[k].mHist_v[0]->GetBinContent(j) 
01197              <<" +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<" )"<<endl;
01198         cout << "##### 4-part v" << j << " = (" 
01199              << histFull[k].mHist_v[1]->GetBinContent(j) 
01200              <<" +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<" )"<<endl;
01201       }
01202       
01203       delete histOfMeanIntegV; 
01204       delete histOfMeanIntegV3;
01205       
01206     } else if (m_M==2) {
01207       
01208       TH1D* histOfMeanIntegV2 = new TH1D(*(histFull[k].mHist_v[0]));
01209       histOfMeanIntegV2->Reset();
01210       
01211       TH1D* histOfMeanIntegV4 = new TH1D(*(histFull[k].mHist_v[1]));
01212       histOfMeanIntegV4->Reset();
01213       
01214       for (int j = 1; j < Flow::nHars+1; j++) {
01215         histOfMeanIntegV2->SetBinContent(j, 1./(meanIntegV2[j-1]*perCent));
01216         histOfMeanIntegV2->SetBinError(j,0.);
01217         histOfMeanIntegV4->SetBinContent(j, -0.5/(meanIntegV4[j-1]*perCent));
01218         histOfMeanIntegV4->SetBinError(j,0.);
01219       }
01220       histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV2);
01221       histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV4);
01222 
01223       //force to be zero if the value is "nan"
01224       for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++){
01225         for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
01226           if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
01227                  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ){
01228             histFull[k].mHist_v[ord]->SetBinContent(j,0.);
01229             histFull[k].mHist_v[ord]->SetBinError(j,0.);
01230           }
01231         }
01232       }
01233 
01234       for (int j = 1; j < Flow::nHars+1; j++) {
01235         cout << "##### 2-part v" << j << " = (" 
01236              << histFull[k].mHist_v[0]->GetBinContent(j) 
01237              <<") +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<endl;
01238         cout << "##### 4-part v" << j << " = (" 
01239              << histFull[k].mHist_v[1]->GetBinContent(j) 
01240              <<") +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<endl;
01241       }
01242       
01243       delete histOfMeanIntegV2; 
01244       delete histOfMeanIntegV4;
01245     }
01246     
01247     for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++)
01248       savedHistNames->AddLast(histFull[k].mHist_v[ord]);
01249     
01250   }
01251 
01252   // GetHistList()->ls();
01253   
01254   // Write most histograms
01255   TFile histFile("flow.cumulant.root", "RECREATE");
01256   for (int k = 0; k < Flow::nSels; k++) {
01257     for (int j = 0; j < Flow::nHars; j++) {
01258       for (int pq  = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) { 
01259         XpqYpqDenomNames->AddLast(histFull[k].histFullHar[j].mCumulG0MixDenom[pq]);
01260       }
01261     }
01262   }
01263 
01264   TVectorD* cumulConstants = new TVectorD(30);
01265   (*cumulConstants)(0)=double(Flow::nHars);
01266   (*cumulConstants)(1)=double(Flow::nSels);
01267   (*cumulConstants)(2)=double(Flow::nSubs);
01268   (*cumulConstants)(3)=double(Flow::nPhiBins);
01269   (*cumulConstants)(4)=double(Flow::nPhiBinsFtpc);
01270   (*cumulConstants)(5)=double(mNEtaBins);
01271   (*cumulConstants)(6)=double(nPtBinsPart);
01272   (*cumulConstants)(7)=double(Flow::nCumulIntegOrders);
01273   (*cumulConstants)(8)=double(Flow::nCumulInteg_qMax);
01274   (*cumulConstants)(9)=double(Flow::nCumulDiffOrders);
01275   (*cumulConstants)(10)=double(Flow::nCumulDiff_qMax);
01276   (*cumulConstants)(11)=double(Flow::nCumulMixHar_pMax);
01277   (*cumulConstants)(12)=double(Flow::nCumulMixHar_qMax);
01278   (*cumulConstants)(13)=r0;
01279   (*cumulConstants)(14)=m_M;
01280   (*cumulConstants)(15)=(strlen(pFlowSelect->PidPart()) != 0) ? 1 : 0;//1,pidflow
01281   (*cumulConstants)(16)=r0Mix;
01282   (*cumulConstants)(17)=double(profScale);
01283 
01284   cumulConstants->Write("CumulConstants",TObject::kOverwrite | TObject::kSingleKey);
01285   savedHistNames->AddLast(cumulConstants);
01286 
01287   TObjString* cumulMethodTag = new TObjString( "cumulNew" );
01288   cumulMethodTag->Write("CumulMethodTag",TObject::kOverwrite | TObject::kSingleKey);
01289   savedHistNames->AddLast(cumulMethodTag);
01290 
01291   savedHistNames->Write();
01292   
01293   histFile.Close();
01294   
01295   // write profile for <G>, the denominator of XpqYpq.
01296   TFile XpqYpqDenomNewFile("denominatorNew.root","RECREATE"); 
01297   XpqYpqDenomNames->Write();
01298   XpqYpqDenomNewFile.Close();
01299   delete XpqYpqDenomNames;
01300     
01301   delete pFlowSelect;
01302     
01303   return StMaker::Finish();
01304 }
01305 
01306 //-----------------------------------------------------------------------
01307 
01308 void StFlowCumulantMaker::SetHistoRanges(Bool_t ftpc_included) {
01309 
01310   if (ftpc_included) {
01311     mEtaMin = Flow::etaMin;
01312     mEtaMax = Flow::etaMax;
01313     mNEtaBins = Flow::nEtaBins;
01314   }
01315   else {
01316     mEtaMin = Flow::etaMinTpcOnly;
01317     mEtaMax = Flow::etaMaxTpcOnly;
01318     mNEtaBins = Flow::nEtaBinsTpcOnly;
01319   }
01320   
01321   return;
01322 }
01323 
01325 //
01326 // $Log: StFlowCumulantMaker.cxx,v $
01327 // Revision 1.22  2006/02/22 19:12:29  posk
01328 // Reduced number of output histograms.
01329 //
01330 // Revision 1.21  2005/10/27 17:53:18  aihong
01331 // changes made so that the cumulant method won't pick up the raw phi weight
01332 //
01333 // Revision 1.20  2004/12/17 15:50:09  aihong
01334 // check in v1{3} code
01335 //
01336 // Revision 1.19  2004/12/09 23:47:08  posk
01337 // Minor changes in code formatting.
01338 // Added hist for TPC primary dca to AnalysisMaker.
01339 //
01340 // Revision 1.18  2004/11/16 21:22:22  aihong
01341 // removed old cumulant method
01342 //
01343 // Revision 1.17  2004/08/24 20:22:39  oldi
01344 // Minor modifications to avoid compiler warnings.
01345 //
01346 // Revision 1.16  2003/09/02 17:58:10  perev
01347 // gcc 3.2 updates + WarnOff
01348 //
01349 // Revision 1.15  2003/07/07 21:58:19  posk
01350 // Made units of momentum GeV/c instead of GeV.
01351 //
01352 // Revision 1.14  2003/03/03 16:24:36  aihong
01353 // blow up 4-part cumulant by 1000 in order to let error bars calculated by ROOT
01354 //
01355 // Revision 1.13  2003/01/10 16:40:36  oldi
01356 // Several changes to comply with FTPC tracks:
01357 // - Switch to include/exclude FTPC tracks introduced.
01358 //   The same switch changes the range of the eta histograms.
01359 // - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
01360 // - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
01361 //   West, FarWest (depending on vertex.z()).
01362 // - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
01363 // - Cut to exclude mu-events with no primary vertex introduced.
01364 //   (This is possible for UPC events and FTPC tracks.)
01365 // - Global DCA cut for FTPC tracks added.
01366 // - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
01367 // - Charge cut for FTPC tracks added.
01368 //
01369 // Revision 1.12  2002/05/19 18:58:00  aihong
01370 // speed up cumulant
01371 //
01372 // Revision 1.11  2002/02/19 14:42:18  jeromel
01373 // Added float.h for Linux 7.2
01374 //
01375 // Revision 1.10  2002/02/02 01:10:12  posk
01376 // Added documentation
01377 //
01378 // Revision 1.9  2002/01/24 20:53:51  aihong
01379 // add histograms for combining results from parallel jobs
01380 //
01381 // Revision 1.8  2002/01/14 23:42:36  posk
01382 // Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
01383 //
01384 // Revision 1.7  2001/12/21 17:01:59  aihong
01385 // minor changes
01386 //
01387 // Revision 1.6  2001/12/18 19:27:27  posk
01388 // "proton" and "antiproton" replaced by "pr+" and "pr-".
01389 //
01390 // Revision 1.5  2001/12/11 22:04:01  posk
01391 // Four sets of phiWgt histograms.
01392 // StFlowMaker StFlowEvent::PhiWeight() changes.
01393 // Cumulant histogram names changed.
01394 //
01395 // Revision 1.4  2001/12/06 01:21:14  jeromel
01396 // Mandatory correction : Extraneous comma removed.
01397 //
01398 // Revision 1.3  2001/11/09 21:14:50  posk
01399 // Switched from CERNLIB to TMath. Using global dca instead of dca.
01400 //
01401 // Revision 1.2  2001/11/08 03:12:24  aihong
01402 // clean up redundant histograms
01403 //
01404 // Revision 1.1  2001/11/02 04:47:42  aihong
01405 // install cumulant maker
01406 //