StRoot: StFlowAnalysisMaker/StFlowDirectCumulantMaker.cxx Source File

00001 
00002 //
00003 // $Id: StFlowDirectCumulantMaker.cxx,v 1.2 2010/06/10 16:33:53 posk Exp $
00004 //
00005 // Authors: Dhevan Gangadharan, UCLA, Dec 2009
00006 //
00008 //
00009 // Description:  Maker to analyze Flow using direct cumulants.
00010 //  It needs only one pass through the data.
00011 //  The outputs are 2D histograms for the needed terms, of events vs centrality and pt.
00012 //  The output file is flow.dirCumulant.root.
00013 //  The macro directCumulants_v2.C then converts these terms to cumulants
00014 //  and then v2(pt) for each centrality and min bias.
00015 //  cent=1 is the most central, cent=0 is the yield weighted average of all the others.
00016 //  The integrated v2 is also calculated by integrating over pt.
00017 //  Some documentation is at http://www.star.bnl.gov/protected/bulkcorr/posk/direct4part/index.html
00018 //
00020 
00021 #include <Stiostream.h>
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #include "StMaker.h"
00025 #include "StFlowDirectCumulantMaker.h"
00026 #include "StFlowMaker/StFlowMaker.h"
00027 #include "StFlowMaker/StFlowEvent.h"
00028 #include "StFlowMaker/StFlowConstants.h"
00029 #include "StFlowMaker/StFlowSelection.h"
00030 #include "PhysicalConstants.h"
00031 #include "SystemOfUnits.h"
00032 #include "TFile.h"
00033 #include "TString.h"
00034 #include "TH2.h"
00035 #include "StMessMgr.h"
00036 #define PR(x) cout << "##### FlowAnalysis: " << (#x) << " = " << (x) << endl;
00037 
00038 ClassImp(StFlowDirectCumulantMaker)
00039 
00040 //-----------------------------------------------------------------------
00041 StFlowDirectCumulantMaker::StFlowDirectCumulantMaker(const Char_t* name): StMaker(name),
00042   MakerName(name) {
00043   pFlowSelect = new StFlowSelection();
00044 }
00045 
00046 StFlowDirectCumulantMaker::StFlowDirectCumulantMaker(const Char_t* name,
00047                                          const StFlowSelection& flowSelect) :
00048   StMaker(name), MakerName(name) {
00049   pFlowSelect = new StFlowSelection(flowSelect); // copy constructor
00050 }
00051 
00052 StFlowDirectCumulantMaker::~StFlowDirectCumulantMaker() {
00053 }
00054 
00055 //-----------------------------------------------------------------------
00056 Int_t StFlowDirectCumulantMaker::Make() {
00057   // do events
00058 
00059   // Get a pointer to StFlowEvent
00060   StFlowMaker* pFlowMaker = NULL;
00061   pFlowMaker = (StFlowMaker*)GetMaker("Flow");
00062   if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
00063   if (pFlowEvent) {
00064     int runID = pFlowEvent->RunID();
00065     //cout << "runID= " << runID << endl;
00066 
00067     if(runID < 8000000) { // not year 7, Mweights = 1
00068       if (Mweights[0][0] != 1.) { // first event
00069         gMessMgr->Info("##### FlowDirCumulant: Setting Multiplicity Weights to 1");
00070         //cout << "Multiplicity Weights:" << endl;
00071         for(int i=0; i<30; i++){
00072           //cout << endl;
00073           for(int j=0; j<1000; j++){
00074             //cout << setprecision(3) << Mweights[i][j] << ", ";        
00075             Mweights[i][j]=1;
00076           }
00077           //cout << endl;
00078         }
00079       }
00080     }
00081 
00082     FillParticleArrays(); // This fills the particle arrays from pFlowTrack
00083     CalculateTerms();
00084     Fill_Histograms(); // Fill Histograms and let things accumulate
00085   } else {
00086     gMessMgr->Info("##### FlowDirCumulant: FlowEvent pointer null");
00087     return kStOK;
00088   }
00089     
00090   if (Debug()) StMaker::PrintInfo();
00091   
00092   return kStOK;
00093 }
00094 
00095 //-----------------------------------------------------------------------
00096 Int_t StFlowDirectCumulantMaker::Init() {
00097   // Book histograms
00098   cout << endl <<"StFlowDirectCumulantMaker::Init()" << endl;
00099 
00100   GetMweights(); // Y7 mutiplicity weights; due to Y7 centrality bias
00101 
00102   TERM1=0;
00103   TERM2=1;
00104   TERM3=2;
00105   TERM4=3;
00106   TERM5=4;
00107   TERM6=5;
00108   TERM7=6;
00109   TERM8=7;
00110   TERM9=8;
00111   TERM10=9;
00112 
00113   DIF=0;
00114   INT=1; 
00115   COS=0;
00116   SIN=1;
00117   
00118   TString *name;
00119   const int ptbins = Flow::PTBINS - 1;
00120   
00121   for(int term=0; term<Flow::TERMS; term++){
00122     for(int type=0; type<Flow::TYPES; type++){
00123       for(int phase=0; phase<Flow::PHASES; phase++){
00124         for(int species=0; species<Flow::SPECIES; species++){
00125                   
00126           name = new TString("Term_");
00127           *name += term+1;
00128 
00129           if(type==0) name->Append("_D");
00130           else name->Append("_I");
00131           
00132           if(phase==0) name->Append("_cos_");
00133           else name->Append("_sin_");
00134 
00135           *name += species;
00136 
00137           Term[term].Type[type].Phase[phase].Species[species].Sum_Singles = new TH2D(name->Data(),"Correlation Sum",Flow::nCents,.5,Flow::nCents+0.5,ptbins,0.05,ptbins/10.+0.05);
00138           Term[term].Type[type].Phase[phase].Species[species].Sum_Singles->GetXaxis()->SetTitle("Centrality bin");
00139           Term[term].Type[type].Phase[phase].Species[species].Sum_Singles->GetYaxis()->SetTitle("p_{t} bin/10");
00140 
00141           name->Append("_Sq");
00142           Term[term].Type[type].Phase[phase].Species[species].Sum_Squares = new TH2D(name->Data(),"Squares Sum",Flow::nCents,.5,Flow::nCents+0.5,ptbins,0.05,ptbins/10.+0.05);
00143           Term[term].Type[type].Phase[phase].Species[species].Sum_Squares->GetXaxis()->SetTitle("Centrality bin");
00144           Term[term].Type[type].Phase[phase].Species[species].Sum_Squares->GetYaxis()->SetTitle("p_{t} bin/10");
00145 
00146         }
00147       }
00148     }
00149   }
00150 
00151   for(int species=0; species<Flow::SPECIES; species++){
00152     name = new TString("Cent_Pt_yields_");
00153     *name += species;
00154     Entries[species].cent_yields = new TH2D(name->Data(),"Pt_spectrum",Flow::nCents,.5,Flow::nCents+0.5,ptbins,0.05,ptbins/10.+0.05);
00155     Entries[species].cent_yields->GetXaxis()->SetTitle("Centrality bin");
00156     Entries[species].cent_yields->GetYaxis()->SetTitle("p_{t} bin");
00157   }
00158   
00159   Event_counter = new TH1D("Event_counter","Event_counter",Flow::nCents,0.5,Flow::nCents+0.5);
00160   Event_counterWeighted = new TH1D("Event_counterWeighted","Event_counterWeighted",Flow::nCents,0.5,Flow::nCents+0.5);
00161   Event_counterWeighted->Sumw2();
00162 
00163   gMessMgr->SetLimit("##### FlowDirCumu", 2);
00164   gMessMgr->Info("##### FlowDirCumu: $Id: StFlowDirectCumulantMaker.cxx,v 1.2 2010/06/10 16:33:53 posk Exp $");
00165 
00166   return StMaker::Init();
00167 }
00168 
00169 //-----------------------------------------------------------------------
00170 void StFlowDirectCumulantMaker::GetMweights() {
00171   // Y7 mutiplicity weights; due to Y7 centrality bias
00172 
00173   TFile* WeightFile = new TFile("M_WeightsY7.root","READ");
00174 
00175   if(WeightFile->IsOpen()) {    
00176     for(int i=0; i<30; i++){ // vertex z
00177       for(int j=0; j<1000; j++){ // ref mult
00178         
00179         Mweights[i][j] = ((TH2D*)WeightFile->Get("Weights"))->GetBinContent(i+1,j+1);
00180         
00181       }
00182     }
00183     WeightFile->Close();
00184   } else {
00185     gMessMgr->Info("##### FlowDirCumulant: No year 7 Multiplicity Weights");
00186     for(int i=0; i<30; i++){
00187       for(int j=0; j<1000; j++){
00188         Mweights[i][j]=1;
00189       }
00190     }
00191   }
00192 }
00193 
00194 //----------------------------------------------------------------------
00195 void StFlowDirectCumulantMaker::FillParticleArrays() {
00196   // Fill _p primary quantities
00197 
00198   grefmult = pFlowEvent->MultEta();
00199 
00200   StThreeVectorF vertex = pFlowEvent->VertexPos();
00201   zbin=0; // z_vertex bin used for multiplicity weights for the Y7 AuAu Centrality bias
00202   for(int i=0; i<30; i++){
00203     if( (vertex[2] > (2*i-30)) && (vertex[2] < (2*(i+1)-30)) ){
00204       zbin=i+1;
00205       break;
00206     }
00207   }
00208 
00209   // Clear previous event quantities
00210   p_count=0;
00211 
00212   for(int i=0; i<Flow::PTBINS; i++) {
00213     pt_bin_counter_p[i]=0;
00214   } 
00215 
00216   // Primary tracks
00217   StFlowTrackCollection* pFlowTracks_primary = pFlowEvent->TrackCollection();
00218   StFlowTrackIterator itr;
00219 
00221   // Primary track arrays
00222   for (itr = pFlowTracks_primary->begin(); itr != pFlowTracks_primary->end(); itr++) {
00223     StFlowTrack* pFlowTrack = *itr;
00224     // cuts were made as StFlowTrack was filled
00225 
00226     phi_p[p_count] = pFlowTrack->Phi();
00227     if(phi_p[p_count] < 0.) phi_p[p_count] += twopi;
00228     id_p[p_count] = pFlowTrack->id();
00229     pt_p[p_count] = pFlowTrack->Pt();
00230 
00231     weights_p[p_count] = 1.0;
00232     if (pFlowEvent->PtWgt()) { // pt wgt
00233       weights_p[p_count] *= (pt_p[p_count] < pFlowEvent->PtWgtSaturation()) ? pt_p[p_count] : 
00234         pFlowEvent->PtWgtSaturation();  // pt weighting going constant
00235     }
00236 
00237     if(pt_p[p_count]      <= .2 ) { ptbin_p[p_count] = 1; goto ptbin; }
00238     else if(pt_p[p_count] <= .4 ) { ptbin_p[p_count] = 3; goto ptbin; } 
00239     else if(pt_p[p_count] <= .6 ) { ptbin_p[p_count] = 5; goto ptbin; } 
00240     else if(pt_p[p_count] <= .8 ) { ptbin_p[p_count] = 7; goto ptbin; }
00241     else if(pt_p[p_count] <= 1.0) { ptbin_p[p_count] = 9; goto ptbin; }
00242     else if(pt_p[p_count] <= 1.2) { ptbin_p[p_count] =11; goto ptbin; }
00243     else if(pt_p[p_count] <= 1.4) { ptbin_p[p_count] =13; goto ptbin; }
00244     else if(pt_p[p_count] <= 1.6) { ptbin_p[p_count] =15; goto ptbin; }
00245     else if(pt_p[p_count] <= 1.8) { ptbin_p[p_count] =17; goto ptbin; }
00246     else if(pt_p[p_count] <= 2.0) { ptbin_p[p_count] =19; goto ptbin; }
00247     else if(pt_p[p_count] <= 2.2) { ptbin_p[p_count] =21; goto ptbin; }
00248     else if(pt_p[p_count] <= 2.4) { ptbin_p[p_count] =23; goto ptbin; }
00249     else if(pt_p[p_count] <= 2.6) { ptbin_p[p_count] =25; goto ptbin; }
00250     else if(pt_p[p_count] <= 2.8) { ptbin_p[p_count] =27; goto ptbin; }
00251     else if(pt_p[p_count] <= 3.2) { ptbin_p[p_count] =30; goto ptbin; }
00252     else if(pt_p[p_count] <= 3.6) { ptbin_p[p_count] =34; goto ptbin; }
00253     else if(pt_p[p_count] <= 4.0) { ptbin_p[p_count] =38; goto ptbin; }
00254     else if(pt_p[p_count] <= 4.4) { ptbin_p[p_count] =42; goto ptbin; }
00255     else if(pt_p[p_count] <= 4.8) { ptbin_p[p_count] =46; goto ptbin; }
00256     else if(pt_p[p_count] <= 5.2) { ptbin_p[p_count] =50; goto ptbin; }
00257     else if(pt_p[p_count] <= 5.6) { ptbin_p[p_count] =54; goto ptbin; }
00258     else if(pt_p[p_count] <= 6.0) { ptbin_p[p_count] =58; goto ptbin; }
00259     else                          { ptbin_p[p_count] =60; }
00260 
00261   ptbin:
00262     pt_bin_counter_p[ptbin_p[p_count]]++;
00263     p_count++;
00264     if(p_count == Flow::MAXMULT) {
00265       gMessMgr->Info("##### FlowDirCumulant: MAXMULT not big enough");
00266       break;
00267     }
00268   }
00269  
00270 }
00271 
00272 //-----------------------------------------------------------------------
00273 void StFlowDirectCumulantMaker::CalculateTerms() {
00274   // Calculate terms
00275 
00276   if(p_count < 4) return; // Not a good event
00277 
00278   ClearTerms(); // Clears all the sum terms
00279   CalculateSumTerms(kTRUE); // Fills the sum terms for the case of integrated flow
00280   CalculateCorrelationsIntegrated(); // Calculates all the integrated correlations: 4-particle, 3-particle, 2-particle, and 1-particle of all needed types
00281 
00282   ClearTerms(); // Clears all the sum terms
00283   CalculateSumTerms(kFALSE);// Fills the sum terms for the case of differential flow
00284   CalculateCorrelationsDifferential(); // Calculates all the differential correlations: 4-particle, 3-particle, 2-particle,and 1-particle of all needed types
00285 
00286 }
00287 
00288 //----------------------------------------------------------------------
00289 void StFlowDirectCumulantMaker::CalculateSumTerms(bool integrated) {
00290   // Calculate component sum terms
00291 
00292   double phi_one   = -1;
00293   int    ptbin_one = -1;
00294   int    id_one    = -1;
00295   double W_one     =  1;
00296   double W_two     =  1;
00297  
00298   // species=0 == charged_hadron,   
00299   // All of the sum terms used below are components of the 4,3,2 particle correlations done with 2 loops. 
00300   // 1-particle correlations are of course still done with 1 loop.
00301 
00302   for(int species=0; species<Flow::SPECIES; species++){
00304     
00305     int max = p_count;
00306     for(int one=0; one<max; one++){
00308       
00309       phi_one   = phi_p[one];
00310       ptbin_one = ptbin_p[one];
00311       id_one    = id_p[one];
00312       if(integrated == kTRUE) W_one = weights_p[one];
00313       else W_one = 1.;
00314        
00315       sum_1_cos_diff[species][ptbin_one] += cos(2*phi_one)*W_one;
00316       sum_1_sin_diff[species][ptbin_one] += sin(2*phi_one)*W_one;
00317       sum_2_cos_diff[species][ptbin_one] += cos(-2*phi_one)*W_one;
00318       sum_2_sin_diff[species][ptbin_one] += sin(-2*phi_one)*W_one;
00319       
00320       sum_1_cos_full[species] += cos(2*phi_one)*W_one;
00321       sum_1_sin_full[species] += sin(2*phi_one)*W_one;
00322       sum_2_cos_full[species] += cos(-2*phi_one)*W_one;
00323       sum_2_sin_full[species] += sin(-2*phi_one)*W_one;      
00324 
00325       for(int two=0; two<p_count; two++){
00327         
00328         if(integrated == kTRUE || species==0) {if(id_p[two] == id_one) continue;}       
00329         W_two = weights_p[two];
00330 
00331         sum_3_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*W_one*W_two;
00332         sum_3_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*W_one*W_two;
00333 
00334         sum_4_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*(phi_one-phi_p[two]))*W_one*W_two*W_one*W_two;
00335         sum_4_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*sin(2*(phi_one-phi_p[two]))*W_one*W_two*W_one*W_two;
00336 
00337         sum_5_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*W_one*W_two*W_two;
00338         sum_5_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*W_one*W_two*W_two;
00339 
00340         sum_6_cos_sin_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*W_one*W_two*W_two;
00341         sum_6_sin_cos_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*W_one*W_two*W_two;
00342 
00343         sum_7_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*cos(2*phi_one)*W_one*W_two*W_one*W_two;
00344 
00345         sum_8_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*cos(2*phi_p[two])*W_one*W_two*W_two*W_two;
00346 
00347         sum_9_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*sin(-2*phi_one)*W_one*W_two*W_one*W_two;
00348 
00349         sum_10_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*sin(-2*phi_p[two])*W_one*W_two*W_two*W_two;
00350 
00351         sum_11_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one+phi_p[two]))*W_one*W_two;
00352         sum_11_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one+phi_p[two]))*W_one*W_two;
00353 
00354       }
00355     }
00356   
00357     int min_bin=0, max_bin=0;
00358     for(int ptbin1=0; ptbin1<Flow::PTBINS; ptbin1++){
00360       if(integrated == kTRUE) {min_bin=0; max_bin=Flow::PTBINS;} 
00361       else {min_bin = ptbin1; max_bin = ptbin1+1;}
00362       
00363       for(int ptbin2=min_bin; ptbin2<max_bin; ptbin2++){
00365         if(integrated == kTRUE) {if(species==0) if(ptbin1==ptbin2) continue;}
00366         
00367         for(int ptbin3=0; ptbin3<Flow::PTBINS; ptbin3++){
00369           if(species==0) if(ptbin1==ptbin3) continue;
00370           
00371           sum_3_cos[species][ptbin1] += sum_3_cos_diff[species][ptbin2][ptbin3]; 
00372           sum_3_sin[species][ptbin1] += sum_3_sin_diff[species][ptbin2][ptbin3];
00373                   
00374           sum_4_cos[species][ptbin1] += sum_4_cos_diff[species][ptbin2][ptbin3];
00375           sum_4_sin[species][ptbin1] += sum_4_sin_diff[species][ptbin2][ptbin3];
00376           
00377           sum_5_cos[species][ptbin1] += sum_5_cos_diff[species][ptbin2][ptbin3];
00378           sum_5_sin[species][ptbin1] += sum_5_sin_diff[species][ptbin2][ptbin3];
00379           
00380           sum_6_cos_sin[species][ptbin1] += sum_6_cos_sin_diff[species][ptbin2][ptbin3];
00381           sum_6_sin_cos[species][ptbin1] += sum_6_sin_cos_diff[species][ptbin2][ptbin3];
00382           
00383           sum_7[species][ptbin1] += sum_7_diff[species][ptbin2][ptbin3];
00384 
00385           sum_8[species][ptbin1] += sum_8_diff[species][ptbin2][ptbin3];
00386 
00387           sum_9[species][ptbin1] += sum_9_diff[species][ptbin2][ptbin3];
00388 
00389           sum_10[species][ptbin1] += sum_10_diff[species][ptbin2][ptbin3];
00390 
00391           sum_11_cos[species][ptbin1] += sum_11_cos_diff[species][ptbin2][ptbin3];
00392           sum_11_sin[species][ptbin1] += sum_11_sin_diff[species][ptbin2][ptbin3];
00393 
00394         }
00395       }
00396    
00397       if(integrated == kTRUE){
00398         sum_3_cos_reference[species][ptbin1] = sum_3_cos[species][ptbin1];
00399         sum_3_sin_reference[species][ptbin1] = sum_3_sin[species][ptbin1];
00400 
00401         sum_11_cos_reference[species][ptbin1] = sum_11_cos[species][ptbin1];
00402         sum_11_sin_reference[species][ptbin1] = sum_11_sin[species][ptbin1];
00403       
00404         if(species==0){
00405           sum_1_cos[species][ptbin1] = sum_1_cos_full[species] - sum_1_cos_diff[species][ptbin1];
00406           sum_1_sin[species][ptbin1] = sum_1_sin_full[species] - sum_1_sin_diff[species][ptbin1];
00407           sum_2_cos[species][ptbin1] = sum_2_cos_full[species] - sum_2_cos_diff[species][ptbin1];
00408           sum_2_sin[species][ptbin1] = sum_2_sin_full[species] - sum_2_sin_diff[species][ptbin1];
00409           
00410           sum_1_cos_reference[species][ptbin1] = sum_1_cos_full[species] - sum_1_cos_diff[species][ptbin1];
00411           sum_1_sin_reference[species][ptbin1] = sum_1_sin_full[species] - sum_1_sin_diff[species][ptbin1];
00412           sum_2_cos_reference[species][ptbin1] = sum_2_cos_full[species] - sum_2_cos_diff[species][ptbin1];
00413           sum_2_sin_reference[species][ptbin1] = sum_2_sin_full[species] - sum_2_sin_diff[species][ptbin1];
00414         }
00415         else {
00416           sum_1_cos[species][ptbin1] = sum_1_cos_full[species];
00417           sum_1_sin[species][ptbin1] = sum_1_sin_full[species];
00418           sum_2_cos[species][ptbin1] = sum_2_cos_full[species];
00419           sum_2_sin[species][ptbin1] = sum_2_sin_full[species];
00420         
00421           sum_1_cos_reference[species][ptbin1] = sum_1_cos_full[species];
00422           sum_1_sin_reference[species][ptbin1] = sum_1_sin_full[species];
00423           sum_2_cos_reference[species][ptbin1] = sum_2_cos_full[species];
00424           sum_2_sin_reference[species][ptbin1] = sum_2_sin_full[species];
00425         }
00426       }
00427       else {
00428         sum_1_cos[species][ptbin1] = sum_1_cos_diff[species][ptbin1];
00429         sum_1_sin[species][ptbin1] = sum_1_sin_diff[species][ptbin1];
00430         sum_2_cos[species][ptbin1] = sum_2_cos_diff[species][ptbin1];
00431         sum_2_sin[species][ptbin1] = sum_2_sin_diff[species][ptbin1];
00432       }
00433 
00434     }
00435 
00436   }
00437     
00438 }
00439 
00440 //----------------------------------------------------------------------
00441 void StFlowDirectCumulantMaker::Combinations(Bool_t integrated, int species, int counter) {
00442   // set the number of combinatinos to ref_combos_4,3,2,1,0
00443 
00444   ref_combos_4 = 0;
00445   ref_combos_3 = 0;
00446   ref_combos_2 = 0;
00447   ref_combos_1 = 0;
00448   ref_combos_0 = 0;
00449 
00450   if(integrated==kTRUE) {
00451     if(species==0) {
00452       ref_combos_4 = (p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1)*(p_count-pt_bin_counter_p[counter]-2);
00453       ref_combos_4 *=(p_count-pt_bin_counter_p[counter]-3);
00454       ref_combos_3 = (p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1)*(p_count-pt_bin_counter_p[counter]-2);
00455       ref_combos_2 = (p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1);
00456       ref_combos_1 = (p_count-pt_bin_counter_p[counter]);
00457       ref_combos_0 = (p_count-pt_bin_counter_p[counter]);
00458     }
00459   }
00460   else {
00461     if(species==0) {
00462       ref_combos_4 = pt_bin_counter_p[counter]*(p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1);
00463       ref_combos_4 *= (p_count-pt_bin_counter_p[counter]-2);
00464       ref_combos_3 = pt_bin_counter_p[counter]*(p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1);
00465       ref_combos_2 = pt_bin_counter_p[counter]*(p_count-pt_bin_counter_p[counter]);
00466       ref_combos_1 = pt_bin_counter_p[counter];
00467       ref_combos_0 = (p_count-pt_bin_counter_p[counter]);
00468     }
00469   }
00470 
00471 }
00472 
00473 //----------------------------------------------------------------------
00474 void StFlowDirectCumulantMaker::CalculateCorrelationsDifferential() {
00475   // Calculate terms for a fixed pt
00476 
00477   // species=0 == charged_hadron
00478   // This part puts the component sum terms together to form the relevant correlations
00479 
00480   for(int species=0; species<Flow::SPECIES; species++){
00481     
00482     for(int ptbin=0; ptbin<Flow::PTBINS; ptbin++){
00483       
00484       Combinations(kFALSE,species,ptbin);
00485 
00486       if(ref_combos_4 == 0) continue;
00487       if(ref_combos_3 == 0) continue;
00488       if(ref_combos_2 == 0) continue;
00489       if(ref_combos_1 == 0) continue;
00490       if(ref_combos_0 == 0) continue;
00491 
00492       Cumulant_Terms[TERM1][DIF][COS][species][ptbin] = sum_3_cos[species][ptbin]*sum_3_cos_reference[species][ptbin]; 
00493       Cumulant_Terms[TERM1][DIF][COS][species][ptbin] -= 2*(sum_5_cos[species][ptbin]*sum_2_cos_reference[species][ptbin] - sum_8[species][ptbin]);
00494       Cumulant_Terms[TERM1][DIF][COS][species][ptbin] += 2*(sum_6_cos_sin[species][ptbin]*sum_2_sin_reference[species][ptbin] - sum_10[species][ptbin]);
00495       Cumulant_Terms[TERM1][DIF][COS][species][ptbin] /= double(ref_combos_4);
00496       
00497       Cumulant_Terms[TERM2][DIF][COS][species][ptbin] = sum_3_cos[species][ptbin]/double(ref_combos_2);
00498       Cumulant_Terms[TERM2][DIF][SIN][species][ptbin] = sum_3_sin[species][ptbin]/double(ref_combos_2);
00499       
00500       Cumulant_Terms[TERM4][DIF][COS][species][ptbin] = sum_1_cos[species][ptbin]/double(ref_combos_1);
00501       Cumulant_Terms[TERM4][DIF][SIN][species][ptbin] = sum_1_sin[species][ptbin]/double(ref_combos_1);
00502 
00503       Cumulant_Terms[TERM7][DIF][COS][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_2_cos_reference[species][ptbin] - sum_5_cos[species][ptbin]));
00504       Cumulant_Terms[TERM7][DIF][COS][species][ptbin] -= (sum_3_sin[species][ptbin]*sum_2_sin_reference[species][ptbin] + sum_5_sin[species][ptbin]);
00505       Cumulant_Terms[TERM7][DIF][COS][species][ptbin] /= double(ref_combos_3);
00506 
00507       Cumulant_Terms[TERM7][DIF][SIN][species][ptbin] = (sum_3_sin[species][ptbin]*sum_2_cos_reference[species][ptbin]) + (sum_3_cos[species][ptbin]*sum_2_sin_reference[species][ptbin]);
00508       Cumulant_Terms[TERM7][DIF][SIN][species][ptbin] -= (sum_6_sin_cos[species][ptbin] - sum_6_cos_sin[species][ptbin]);
00509       Cumulant_Terms[TERM7][DIF][SIN][species][ptbin] /= double(ref_combos_3);
00510 
00511       Cumulant_Terms[TERM9][DIF][COS][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_1_cos_reference[species][ptbin] - sum_5_cos[species][ptbin]));
00512       Cumulant_Terms[TERM9][DIF][COS][species][ptbin] -= (sum_3_sin[species][ptbin]*sum_1_sin_reference[species][ptbin] - sum_5_sin[species][ptbin]);
00513       Cumulant_Terms[TERM9][DIF][COS][species][ptbin] /= double(ref_combos_3);      
00514 
00515       Cumulant_Terms[TERM9][DIF][SIN][species][ptbin] = (sum_3_sin[species][ptbin]*sum_1_cos_reference[species][ptbin]) + (sum_3_cos[species][ptbin]*sum_1_sin_reference[species][ptbin]);
00516       Cumulant_Terms[TERM9][DIF][SIN][species][ptbin] -= (sum_6_sin_cos[species][ptbin] + sum_6_cos_sin[species][ptbin]);
00517       Cumulant_Terms[TERM9][DIF][SIN][species][ptbin] /= double(ref_combos_3);
00518 
00519       Cumulant_Terms[TERM10][DIF][COS][species][ptbin] = sum_11_cos[species][ptbin]/double(ref_combos_2);
00520       Cumulant_Terms[TERM10][DIF][SIN][species][ptbin] = sum_11_sin[species][ptbin]/double(ref_combos_2);
00521        
00522     }
00523   }
00524 }
00525 
00526 //----------------------------------------------------------------------
00527 void StFlowDirectCumulantMaker::CalculateCorrelationsIntegrated() {
00528   // Calculate terms missing a pt to avoid autocorrelations
00529     
00530   // species=0 == charged_hadron
00531   // Similar to CalculateCorrelationsDifferential but for the integrated flow case
00532 
00533   for(int species=0; species<Flow::SPECIES; species++){
00534 
00535     for(int ptbin=0; ptbin<Flow::PTBINS; ptbin++){
00536 
00537       Combinations(kTRUE,species,ptbin);
00538 
00539       if(ref_combos_4 == 0) continue;
00540       if(ref_combos_3 == 0) continue;
00541       if(ref_combos_2 == 0) continue;
00542       if(ref_combos_1 == 0) continue;
00543       if(ref_combos_0 == 0) continue;
00544       
00545       Cumulant_Terms[TERM1][INT][COS][species][ptbin] = sum_3_cos[species][ptbin]*sum_3_cos[species][ptbin] - 2*sum_4_cos[species][ptbin];
00546       Cumulant_Terms[TERM1][INT][COS][species][ptbin] -= 4*((sum_5_cos[species][ptbin]*sum_2_cos[species][ptbin] - sum_8[species][ptbin] - sum_7[species][ptbin]) - (sum_6_cos_sin[species][ptbin]*sum_2_sin[species][ptbin] - sum_10[species][ptbin] - sum_9[species][ptbin]));
00547       Cumulant_Terms[TERM1][INT][COS][species][ptbin] /= double(ref_combos_4);
00548    
00549       Cumulant_Terms[TERM2][INT][COS][species][ptbin] = sum_3_cos[species][ptbin]/double(ref_combos_2);
00550       Cumulant_Terms[TERM2][INT][SIN][species][ptbin] = sum_3_sin[species][ptbin]/double(ref_combos_2);
00551       Cumulant_Terms[TERM3][INT][COS][species][ptbin] = sum_3_cos[species][ptbin]/double(ref_combos_2);
00552       Cumulant_Terms[TERM3][INT][SIN][species][ptbin] = sum_3_sin[species][ptbin]/double(ref_combos_2);
00553 
00554       Cumulant_Terms[TERM4][INT][COS][species][ptbin] = sum_1_cos[species][ptbin]/double(ref_combos_1);
00555       Cumulant_Terms[TERM4][INT][SIN][species][ptbin] = sum_1_sin[species][ptbin]/double(ref_combos_1);
00556 
00557       Cumulant_Terms[TERM5][INT][COS][species][ptbin] = (sum_3_cos[species][ptbin]*sum_2_cos[species][ptbin] - sum_3_sin[species][ptbin]*sum_2_sin[species][ptbin]);
00558       Cumulant_Terms[TERM5][INT][COS][species][ptbin] -= 2*(sum_5_cos[species][ptbin]);
00559       Cumulant_Terms[TERM5][INT][COS][species][ptbin] /= double(ref_combos_3);
00560       Cumulant_Terms[TERM5][INT][SIN][species][ptbin] = (sum_3_cos[species][ptbin]*sum_2_sin[species][ptbin] + sum_3_sin[species][ptbin]*sum_2_cos[species][ptbin]);
00561 
00562       Cumulant_Terms[TERM5][INT][SIN][species][ptbin] += 2*(sum_6_cos_sin[species][ptbin]);
00563       Cumulant_Terms[TERM5][INT][SIN][species][ptbin] /= double(ref_combos_3);
00564 
00565       Cumulant_Terms[TERM6][INT][COS][species][ptbin] = sum_1_cos[species][ptbin]/double(ref_combos_1);
00566       Cumulant_Terms[TERM6][INT][SIN][species][ptbin] = sum_1_sin[species][ptbin]/double(ref_combos_1);
00567 
00568       Cumulant_Terms[TERM7][INT][COS][species][ptbin] = Cumulant_Terms[TERM5][INT][COS][species][ptbin];
00569       Cumulant_Terms[TERM7][INT][SIN][species][ptbin] = Cumulant_Terms[TERM5][INT][SIN][species][ptbin];
00570   
00571       Cumulant_Terms[TERM8][INT][COS][species][ptbin] = sum_2_cos[species][ptbin]/double(ref_combos_1);
00572       Cumulant_Terms[TERM8][INT][SIN][species][ptbin] = sum_2_sin[species][ptbin]/double(ref_combos_1);
00573 
00574       Cumulant_Terms[TERM9][INT][COS][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_1_cos[species][ptbin] - 2*sum_5_cos[species][ptbin]));
00575       Cumulant_Terms[TERM9][INT][COS][species][ptbin] -= (sum_3_sin[species][ptbin]*sum_1_sin[species][ptbin]);
00576       Cumulant_Terms[TERM9][INT][COS][species][ptbin] /= double(ref_combos_3);
00577       Cumulant_Terms[TERM9][INT][SIN][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_1_sin[species][ptbin] - 2*sum_6_cos_sin[species][ptbin]));
00578       Cumulant_Terms[TERM9][INT][SIN][species][ptbin] += (sum_3_sin[species][ptbin]*sum_1_cos[species][ptbin]);
00579       Cumulant_Terms[TERM9][INT][SIN][species][ptbin] /= double(ref_combos_3);
00580 
00581       Cumulant_Terms[TERM10][INT][COS][species][ptbin] = sum_11_cos[species][ptbin]/double(ref_combos_2);
00582       Cumulant_Terms[TERM10][INT][SIN][species][ptbin] = sum_11_sin[species][ptbin]/double(ref_combos_2);
00583       
00584     }
00585   }
00586 
00587 }
00588 
00589 //----------------------------------------------------------------------
00590 void StFlowDirectCumulantMaker::Fill_Histograms() {
00591   // Fill term histograms
00592   
00593   Event_counter->Fill(pFlowEvent->Centrality());
00594   if (Mweights[zbin-1][grefmult-1] < 12.2) {
00595     Event_counterWeighted->Fill(pFlowEvent->Centrality(), Mweights[zbin-1][grefmult-1]);
00596   }
00597 
00598   for(int ptbin=0; ptbin<Flow::PTBINS; ptbin++){
00599     for(int species=0; species<Flow::SPECIES; species++){
00600             
00601       double the_pt = double(double(ptbin)/10.);  // GeV = bin/10.
00602       if(ptbin == Flow::PTBINS-1) Combinations(kTRUE,species,ptbin);
00603       else Combinations(kFALSE,species,ptbin);
00604 
00605       if(ref_combos_4 == 0) continue;
00606       if(ref_combos_3 == 0) continue;
00607       if(ref_combos_2 == 0) continue;
00608       if(ref_combos_1 == 0) continue;
00609       if(ref_combos_0 == 0) continue;
00610 
00611       if(Mweights[zbin-1][grefmult-1] > 12.2) continue;
00612       double PF = ref_combos_1*Mweights[zbin-1][grefmult-1]; // prefactor for Y7 centrality bias correction
00613 
00614       for(int phase=0; phase<Flow::PHASES; phase++){
00615         for(int type=0; type<Flow::TYPES; type++){
00616           for(int term=0; term<Flow::TERMS; term++){
00617             if(fabs(Cumulant_Terms[term][type][phase][species][ptbin]) < 1.) {
00618               Term[term].Type[type].Phase[phase].Species[species].Sum_Singles->Fill(pFlowEvent->Centrality(),the_pt, PF*Cumulant_Terms[term][type][phase][species][ptbin]);
00619             
00620               Term[term].Type[type].Phase[phase].Species[species].Sum_Squares->Fill(pFlowEvent->Centrality(),the_pt, PF*pow(Cumulant_Terms[term][type][phase][species][ptbin],2));
00621             }
00622             
00623             
00624           }
00625         }
00626       }
00627       Entries[species].cent_yields->Fill(pFlowEvent->Centrality(), the_pt, ref_combos_1*Mweights[zbin-1][grefmult-1]);
00628     }
00629   }
00630   
00631 }
00632 
00633 //----------------------------------------------------------------------
00634 void StFlowDirectCumulantMaker::ClearTerms() {
00635   // All the component sum terms need to be cleared twice per event: 1 for the Differential calculation, 1 for the Integrated calculation.
00636   
00637   for(int species=0; species<Flow::SPECIES; species++){
00638 
00639     sum_1_cos_full[species]=0;
00640     sum_1_sin_full[species]=0;
00641     sum_2_cos_full[species]=0;
00642     sum_2_sin_full[species]=0;
00643 
00644     for(int ptbin1=0; ptbin1<Flow::PTBINS; ptbin1++){
00645 
00646       for(int ptbin2=0; ptbin2<Flow::PTBINS; ptbin2++){
00647         
00648         sum_1_cos_diff[species][ptbin1]=0, sum_1_sin_diff[species][ptbin1]=0;
00649         sum_2_cos_diff[species][ptbin1]=0, sum_2_sin_diff[species][ptbin1]=0;
00650         sum_3_cos_diff[species][ptbin1][ptbin2]=0, sum_3_sin_diff[species][ptbin1][ptbin2]=0;
00651         sum_4_cos_diff[species][ptbin1][ptbin2]=0, sum_4_sin_diff[species][ptbin1][ptbin2]=0;
00652         sum_5_cos_diff[species][ptbin1][ptbin2]=0, sum_5_sin_diff[species][ptbin1][ptbin2]=0;
00653         sum_6_cos_sin_diff[species][ptbin1][ptbin2]=0, sum_6_sin_cos_diff[species][ptbin1][ptbin2]=0;
00654         sum_7_diff[species][ptbin1][ptbin2]=0, sum_8_diff[species][ptbin1][ptbin2]=0, sum_9_diff[species][ptbin1][ptbin2]=0, sum_10_diff[species][ptbin1][ptbin2]=0;
00655         sum_11_cos_diff[species][ptbin1][ptbin2]=0, sum_11_sin_diff[species][ptbin1][ptbin2]=0;
00656 
00657         sum_1_cos[species][ptbin1]=0, sum_1_sin[species][ptbin1]=0;
00658         sum_2_cos[species][ptbin1]=0, sum_2_sin[species][ptbin1]=0;
00659         sum_3_cos[species][ptbin1]=0, sum_3_sin[species][ptbin1]=0;
00660         sum_4_cos[species][ptbin1]=0, sum_4_sin[species][ptbin1]=0;
00661         sum_5_cos[species][ptbin1]=0, sum_5_sin[species][ptbin1]=0;
00662         sum_6_cos_sin[species][ptbin1]=0, sum_6_sin_cos[species][ptbin1]=0;
00663         sum_7[species][ptbin1]=0, sum_8[species][ptbin1]=0, sum_9[species][ptbin1]=0, sum_10[species][ptbin1]=0;
00664         sum_11_cos[species][ptbin1]=0, sum_11_sin[species][ptbin1]=0;
00665 
00666       }
00667       
00668     }
00669   }
00670   
00671 }
00672 
00673 //-----------------------------------------------------------------------
00674 Int_t StFlowDirectCumulantMaker::Finish() {
00675   // Wtite out file of term histograms
00676 
00677   cout << endl << "##### Direct Cumulant Maker:" << endl;
00678      
00679   // Write all terms
00680   const char *flowname = "flow.dirCumulant.root";
00681   TFile histFile(flowname, "RECREATE");
00682   
00683   GetHistList()->Write();
00684   //GetHistList()->ls();
00685   
00686   histFile.Close();
00687   
00688   delete pFlowSelect;
00689   
00690   return StMaker::Finish();
00691 }
00692 
00693 //-----------------------------------------------------------------------