StRoot: StFlowAnalysisMaker/directCumulants

00001 
00002 //
00003 // $Id: directCumulants_v2.C,v 1.3 2010/06/10 16:33:57 posk Exp $
00004 //
00005 // Authors: Dhevan Gangadharan, UCLA, Jan 2010
00006 //
00008 //
00009 // Description:  Macro to sum the Terms from the output of StFlowDirectCumulantMaker,
00010 //    plot v_2{2] and v_2{4}, and output the graphs to flow.diCrumulant.graphs.root.
00011 //    You can plot them again with plotDirectCumulants_v2pt.C .
00012 //    cent=1 is the most central, cent=0 is the yield weighted average of all the others.
00013 //
00015 
00016 #include <math.h>
00017 #include <stdio.h>
00018 #include <stdlib.h>
00019 #include <Riostream.h>
00020 
00021 #include "TFile.h"
00022 #include "TString.h"
00023 #include "TH1.h"
00024 #include "TH2.h"
00025 #include "TGraphErrors.h"
00026 #include "TArray.h"
00027 #include "TLegend.h"
00028 #include "TStyle.h"
00029 
00030 void directCumulants_v2(){
00031 
00032   TCanvas* cOld = (TCanvas*)gROOT->GetListOfCanvases(); // delete old canvas
00033   if (cOld) cOld->Delete();
00034 
00035   bool ptGraphs = kFALSE;
00036 
00037   //TFile *myfile = new TFile("Y4_test4.root","READ");
00038   TFile *myfile = new TFile("flow.dirCumulant.root","READ");
00039   //TFile *myfile = new TFile("~dhevang/test.root","READ");
00040   const int CENTBINS= 9;
00041   float x[CENTBINS] = {2.5, 7.5, 15., 25., 35., 45., 55., 65., 75.};
00042   
00043   gROOT->SetStyle("Plain");                              // set style
00044   //gStyle->SetOptDate(0);
00045   gStyle->SetOptStat(0);
00046   
00047   const int TERMS=    10; // terms needed to calculate the cumulants
00048   const int TYPES=     2; // differential or integrated
00049   const int PHASES=    2; // sin or cos
00050   const int SPECIES=   1; // charged particles = 0
00051   const int PTBINS=   62; // max pt bins
00052   const int PTINTBINS=19; // max pt bins for pt integration
00053   int ptOverflow = PTBINS-1;
00054   int ptCu = PTBINS-2;    // integrated cumulants
00055   
00056   double yields[SPECIES][PTBINS]={{0}};
00057   double cent_yields[SPECIES][CENTBINS][PTBINS]={{{0}}};
00058   
00059   // centrality 0 is min bias
00060   double v2_2[SPECIES][CENTBINS+1][PTBINS]={{{0.}}}; // differential
00061   double v2_2_e[SPECIES][CENTBINS+1][PTBINS]={{{0.}}};
00062   double v2_4[SPECIES][CENTBINS+1][PTBINS]={{{0.}}};
00063   double v2_4_e[SPECIES][CENTBINS+1][PTBINS]={{{0.}}};
00064   double v2IntPt_2[SPECIES][CENTBINS+1]={{0.}}; // integral
00065   double v2IntPt_2_e[SPECIES][CENTBINS+1]={{0.}};
00066   double v2IntPt_4[SPECIES][CENTBINS+1]={{0.}};
00067   double v2IntPt_4_e[SPECIES][CENTBINS+1]={{0.}};
00068   double used_cent_yields2[SPECIES][CENTBINS+1]={{0.}};
00069   double used_cent_yields4[SPECIES][CENTBINS+1]={{0.}};
00070   double usedInt_cent_yields2[SPECIES][CENTBINS+1]={{0.}};
00071   double usedInt_cent_yields4[SPECIES][CENTBINS+1]={{0.}};
00072  
00073   TString *name;
00074   TString *name1;
00075   TString *name2;
00076   TString *name3;
00077   TString *name4;
00078   
00079   double Terms[TERMS][TYPES][PHASES][SPECIES][CENTBINS][PTBINS]={{{{{{0}}}}}};
00080   double Terms_e[TERMS][TYPES][PHASES][SPECIES][CENTBINS][PTBINS]={{{{{{0}}}}}};
00081   
00083   // This part takes all the info out of the input file and puts it into the "Terms" arrays
00084   
00085   for(int species=0; species<SPECIES; species++){
00086 
00087     name = new TString("Cent_Pt_yields_");
00088     *name += species;
00089 
00090     for(int cent=0; cent<CENTBINS; cent++){
00091       
00092       for(int ptbin=0; ptbin<PTBINS; ptbin++){
00093         
00094         cent_yields[species][cent][ptbin] = ((TH2D*)myfile->Get(name->Data()))->GetBinContent(cent+1,ptbin+1);
00095         yields[species][ptbin] += cent_yields[species][cent][ptbin];
00096         
00097       }
00098     }
00099   }
00100   
00101   for(int term=0; term<TERMS; term++){
00102 
00103     name1 = new TString("Term_");
00104     *name1 += term+1;
00105 
00106     for(int type=0; type<TYPES; type++){
00107       
00108       name2 = new TString();
00109       if(type==0) name2->Append("_D_"); // differential
00110       else name2->Append("_I_"); // integral
00111 
00112       for(int phase=0; phase<PHASES; phase++){
00113 
00114         name3 = new TString();
00115         if(phase==0) name3->Append("cos_");
00116         else name3->Append("sin_");
00117         for(int species=0; species<SPECIES; species++){
00118 
00119           name4 = new TString();
00120           *name4 += species;
00121 
00122           for(int cent=0; cent<CENTBINS; cent++){
00123             
00124             for(int ptbin=0; ptbin<PTBINS; ptbin++){
00125               
00126               if(cent_yields[species][cent][ptbin] < 1) continue;
00127 
00128               name = new TString();
00129               name->Append(name1->Data());
00130               name->Append(name2->Data());
00131               name->Append(name3->Data());
00132               name->Append(name4->Data());
00133 
00134               Terms[term][type][phase][species][cent][ptbin] = ((TH2D*)myfile->Get(name->Data()))->GetBinContent(cent+1,ptbin+1)/cent_yields[species][cent][ptbin];
00135 
00136               name->Append("_Sq");            
00137               Terms_e[term][type][phase][species][cent][ptbin] = ((TH2D*)myfile->Get(name->Data()))->GetBinContent(cent+1,ptbin+1)/cent_yields[species][cent][ptbin];
00138               Terms_e[term][type][phase][species][cent][ptbin] -= pow(Terms[term][type][phase][species][cent][ptbin],2);
00139               Terms_e[term][type][phase][species][cent][ptbin] /= cent_yields[species][cent][ptbin];
00140               Terms_e[term][type][phase][species][cent][ptbin] = sqrt(Terms_e[term][type][phase][species][cent][ptbin]);
00141               
00142             }
00143           }
00144         }
00145       }
00146     }
00147   }
00148 
00149 
00151   int COS=0, SIN=1, I=1;
00152 
00153   double Cumulant_term[11][TYPES][SPECIES][CENTBINS][PTBINS]={{{{{0}}}}};
00154   double Cumulant_term_e[11][TYPES][SPECIES][CENTBINS][PTBINS]={{{{{0}}}}};
00155   double Cumulant_4[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
00156   double Cumulant_4_e[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
00157   double Cumulant_2[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
00158   double Cumulant_2_e[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
00159 
00160   // This part reconstructs all the Cumulant terms from "Terms"
00161   for(int type=0; type<TYPES; type++){
00162     
00163     for(int species=0; species<SPECIES; species++){
00164 
00165       for(int cent=0; cent<CENTBINS; cent++){
00166         
00167         for(int ptbin=0; ptbin<PTBINS; ptbin++){
00168                   
00169           Cumulant_term[0][type][species][cent][ptbin] = Terms[0][type][COS][species][cent][ptbin];
00170           Cumulant_term_e[0][type][species][cent][ptbin] = Terms_e[0][type][COS][species][cent][ptbin];
00171 
00172           Cumulant_term[1][type][species][cent][ptbin] = Terms[1][type][COS][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin];
00173           Cumulant_term[1][type][species][cent][ptbin] -= Terms[1][type][SIN][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin];
00174           Cumulant_term_e[1][type][species][cent][ptbin] = pow(Terms_e[1][type][COS][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin],2);
00175           Cumulant_term_e[1][type][species][cent][ptbin] += pow(Terms[1][type][COS][species][cent][ptbin]*Terms_e[2][I][COS][species][cent][ptbin],2);
00176           Cumulant_term_e[1][type][species][cent][ptbin] += pow(Terms_e[1][type][SIN][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin],2);
00177           Cumulant_term_e[1][type][species][cent][ptbin] += pow(Terms[1][type][SIN][species][cent][ptbin]*Terms_e[2][I][SIN][species][cent][ptbin],2);
00178           Cumulant_term_e[1][type][species][cent][ptbin] = sqrt(Cumulant_term_e[1][type][species][cent][ptbin]);
00179 
00180           Cumulant_term[2][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[4][I][COS][species][cent][ptbin];
00181           Cumulant_term[2][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[4][I][SIN][species][cent][ptbin];
00182           
00183           Cumulant_term[3][type][species][cent][ptbin] = Terms[5][I][COS][species][cent][ptbin]*Terms[6][type][COS][species][cent][ptbin];
00184           Cumulant_term[3][type][species][cent][ptbin] -= Terms[5][I][SIN][species][cent][ptbin]*Terms[6][type][SIN][species][cent][ptbin];
00185 
00186           Cumulant_term[4][type][species][cent][ptbin] = Terms[7][I][COS][species][cent][ptbin]*Terms[8][type][COS][species][cent][ptbin];
00187           Cumulant_term[4][type][species][cent][ptbin] -= Terms[7][I][SIN][species][cent][ptbin]*Terms[8][type][SIN][species][cent][ptbin];
00188 
00189           Cumulant_term[5][type][species][cent][ptbin] = Terms[9][type][COS][species][cent][ptbin]*Terms[9][I][COS][species][cent][ptbin];
00190           Cumulant_term[5][type][species][cent][ptbin] -= Terms[9][type][SIN][species][cent][ptbin]*Terms[9][I][SIN][species][cent][ptbin];
00191 
00192           Cumulant_term[6][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin];
00193           Cumulant_term[6][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin];
00194           Cumulant_term[6][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin];
00195           Cumulant_term[6][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin];
00196 
00197           Cumulant_term_e[6][type][species][cent][ptbin] = pow(Terms[3][type][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms_e[2][I][COS][species][cent][ptbin],2);
00198           Cumulant_term_e[6][type][species][cent][ptbin] = sqrt(Cumulant_term_e[6][I][species][cent][ptbin]);
00199 
00200 
00201           Cumulant_term[7][type][species][cent][ptbin] = Terms[3][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[1][type][COS][species][cent][ptbin];
00202           Cumulant_term[7][type][species][cent][ptbin] -= Terms[3][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[1][type][COS][species][cent][ptbin];
00203           Cumulant_term[7][type][species][cent][ptbin] -= Terms[3][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[1][type][SIN][species][cent][ptbin];
00204           Cumulant_term[7][type][species][cent][ptbin] -= Terms[3][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[1][type][SIN][species][cent][ptbin];
00205 
00206           Cumulant_term_e[7][type][species][cent][ptbin] = pow(Terms[3][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms_e[1][type][COS][species][cent][ptbin],2);
00207           Cumulant_term_e[7][type][species][cent][ptbin] = sqrt(Cumulant_term_e[7][type][species][cent][ptbin]);
00208 
00209 
00210           Cumulant_term[8][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[9][I][COS][species][cent][ptbin];
00211           Cumulant_term[8][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[9][I][COS][species][cent][ptbin];
00212           Cumulant_term[8][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[9][I][SIN][species][cent][ptbin];
00213           Cumulant_term[8][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[9][I][SIN][species][cent][ptbin];
00214           
00215           Cumulant_term[9][type][species][cent][ptbin] = Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[9][type][COS][species][cent][ptbin];
00216           Cumulant_term[9][type][species][cent][ptbin] -= Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[9][type][COS][species][cent][ptbin];
00217           Cumulant_term[9][type][species][cent][ptbin] -= Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[9][type][SIN][species][cent][ptbin];
00218           Cumulant_term[9][type][species][cent][ptbin] -= Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[9][type][SIN][species][cent][ptbin];
00219 
00220           Cumulant_term[10][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
00221           Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
00222           Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
00223           Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
00224           Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
00225           Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
00226           Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
00227           Cumulant_term[10][type][species][cent][ptbin] += Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
00228 
00229 
00230           Cumulant_4[type][species][cent][ptbin] = Cumulant_term[0][type][species][cent][ptbin] -2*Cumulant_term[1][type][species][cent][ptbin];
00231           Cumulant_4[type][species][cent][ptbin] -= Cumulant_term[2][type][species][cent][ptbin] + Cumulant_term[3][type][species][cent][ptbin];
00232           Cumulant_4[type][species][cent][ptbin] -= 2*Cumulant_term[4][type][species][cent][ptbin];
00233           Cumulant_4[type][species][cent][ptbin] -= Cumulant_term[5][type][species][cent][ptbin];
00234           Cumulant_4[type][species][cent][ptbin] += 4*Cumulant_term[6][type][species][cent][ptbin];
00235           Cumulant_4[type][species][cent][ptbin] += 4*Cumulant_term[7][type][species][cent][ptbin];
00236           Cumulant_4[type][species][cent][ptbin] += 2*Cumulant_term[8][type][species][cent][ptbin];
00237           Cumulant_4[type][species][cent][ptbin] += 2*Cumulant_term[9][type][species][cent][ptbin];
00238           Cumulant_4[type][species][cent][ptbin] -= 6*Cumulant_term[10][type][species][cent][ptbin];
00239           
00241 //        if (ptbin==ptCu && type==1) { // integrated cumulants
00242 //          cout << endl << "centrality= " << cent+1 << endl;
00243 //          double noAcc = Cumulant_term[0][type][species][cent][ptbin] -2*Cumulant_term[1][type][species][cent][ptbin];
00244 //          cout << "noAcc: " << setprecision(4) << noAcc << endl;
00245 //          cout << " 2: " << -Cumulant_term[2][type][species][cent][ptbin] << endl;
00246 //          cout << " 3: " << -Cumulant_term[3][type][species][cent][ptbin] << endl;
00247 //          cout << " 4: " << -2*Cumulant_term[4][type][species][cent][ptbin] << endl;
00248 //          cout << " 5: " << -Cumulant_term[5][type][species][cent][ptbin] << endl;
00249 //          cout << " 6: " <<  4*Cumulant_term[6][type][species][cent][ptbin] << endl;
00250 //          cout << " 7: " <<  4*Cumulant_term[7][type][species][cent][ptbin] << endl;
00251 //          cout << " 8: " <<  2*Cumulant_term[8][type][species][cent][ptbin] << endl;
00252 //          cout << " 9: " <<  2*Cumulant_term[9][type][species][cent][ptbin] << endl;
00253 //          cout << "10: " << -6*Cumulant_term[10][type][species][cent][ptbin] << endl;
00254 //          cout << "Acc/noAcc: " << Cumulant_4[type][species][cent][ptbin] / noAcc << endl;
00255 //          if(Cumulant_4[1][0][cent][ptCu]<0 && noAcc<0) {
00256 //            cout << "v2NoAcc: " << pow(-noAcc,0.25) << endl;
00257 //            cout << "v2Acc: " << pow(-Cumulant_4[1][0][cent][ptCu],0.25) << endl;
00258 //            cout << "v2Acc/noAcc: " << pow(-Cumulant_4[1][0][cent][ptCu],0.25) / pow(-noAcc,0.25) << endl;
00259 //          }
00260 //        }
00262           
00263           Cumulant_4_e[type][species][cent][ptbin] = Cumulant_term_e[0][type][species][cent][ptbin]**2;
00264           Cumulant_4_e[type][species][cent][ptbin] += 2*(Cumulant_term_e[1][type][species][cent][ptbin])**2;
00265           Cumulant_4_e[type][species][cent][ptbin] += 4*(Cumulant_term_e[6][type][species][cent][ptbin])**2;
00266           Cumulant_4_e[type][species][cent][ptbin] += 4*(Cumulant_term_e[7][type][species][cent][ptbin])**2;
00267           Cumulant_4_e[type][species][cent][ptbin] = sqrt(Cumulant_4_e[type][species][cent][ptbin]);
00268           
00269 
00270           Cumulant_2[type][species][cent][ptbin] = Terms[1][type][COS][species][cent][ptbin];
00271           Cumulant_2[type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[7][type][COS][species][cent][ptbin];
00272           
00273           Cumulant_2_e[type][species][cent][ptbin] = Terms_e[1][type][COS][species][cent][ptbin]**2;
00274           Cumulant_2_e[type][species][cent][ptbin] += (Terms_e[3][type][COS][species][cent][ptbin]*Terms[7][type][COS][species][cent][ptbin])**2;
00275           Cumulant_2_e[type][species][cent][ptbin] += (Terms[3][type][COS][species][cent][ptbin]*Terms_e[7][type][COS][species][cent][ptbin])**2;
00276           Cumulant_2_e[type][species][cent][ptbin] = sqrt(Cumulant_2_e[type][species][cent][ptbin]);
00277 
00278          
00279         }
00280       }
00281     }
00282   }
00283 
00285   // each pt bin is calculated for all centralities, then the next pt bin
00286   int    ptcount[SPECIES]={0}; // index of pt bins
00287   double used_pt_yields4[CENTBINS+1] = {0.};
00288   double used_pt_yields2[CENTBINS+1] = {0.};
00289   double v2IntCu_2[CENTBINS] = {0.};
00290   double v2IntCu_2Err[CENTBINS] = {0.};
00291   double v2IntCu_4[CENTBINS] = {0.};
00292   double v2IntCu_4Err[CENTBINS] = {0.};
00293   
00294   // This part reconstructs the elliptic flow from the Cumulant terms
00295   for(int species=0; species<SPECIES; species++){
00296     
00297     cout << endl << "cent: \t v2IntCu_2 +/- err,\t v2IntCu_4 +/- err" << endl;
00298 
00299     for(int ptbin=0; ptbin<PTBINS; ptbin++){
00300       
00301       // defaults for min bias set to an unphysical result
00302       v2_2[species][0][ptcount[species]]   = 10.;
00303       v2_2_e[species][0][ptcount[species]] = 10.;
00304       v2_4[species][0][ptcount[species]]   = 10.;
00305       v2_4_e[species][0][ptcount[species]] = 10.;
00306       used_pt_yields4[0] = 0.; // reset for each pt and species
00307       used_pt_yields2[0] = 0.;
00308       
00309       for(int cent=0; cent<CENTBINS; cent++){
00310         if (ptbin==ptCu) { // integrated cumulants
00311           if(Cumulant_2[1][0][cent][ptCu]>0) {
00312             v2IntCu_2[cent]    = pow(Cumulant_2[1][0][cent][ptCu],0.5);
00313             v2IntCu_2Err[cent] = pow((.5/(Cumulant_2[1][0][cent][ptCu])),0.5)*Cumulant_2_e[1][0][cent][ptCu];
00314           }
00315           if(Cumulant_4[1][0][cent][ptCu]<0) {
00316             v2IntCu_4[cent]    = pow(-Cumulant_4[1][0][cent][ptCu],0.25);
00317             v2IntCu_4Err[cent] = pow((.25/(-Cumulant_4[1][0][cent][ptCu])),0.75)*Cumulant_4_e[1][0][cent][ptCu];
00318           }
00319           cout << cent+1 << ": \t" << setprecision(3) << v2IntCu_2[cent]*100. << 
00320             " \t+/- " << v2IntCu_2Err[cent]*100. << ", \t" << v2IntCu_4[cent]*100. <<
00321             " \t+/- " << v2IntCu_4Err[cent]*100. << endl;
00322         }
00323 
00324         // defaults for each cent set to an unphysical result
00325         v2_2[species][cent+1][ptcount[species]]   = 10.;
00326         v2_2_e[species][cent+1][ptcount[species]] = 10.;
00327         v2_4[species][cent+1][ptcount[species]]   = 10.;
00328         v2_4_e[species][cent+1][ptcount[species]] = 10.;
00329         used_pt_yields4[cent+1] = 0.;
00330         used_pt_yields2[cent+1] = 0.;
00331 
00332         // 4-particle cumulant v2
00333         if(Cumulant_4[1][species][cent][ptbin] < 0 && Cumulant_4[0][species][cent][ptbin] < 0) { // both must be negative
00334           v2_4[species][cent+1][ptcount[species]] += -cent_yields[species][cent][ptbin]*Cumulant_4[0][species][cent][ptbin]/(pow(-Cumulant_4[1][species][cent][ptbin],.75));
00335           v2_4[species][0][ptcount[species]] += v2_4[species][cent+1][ptcount[species]];
00336           if(ptbin<PTINTBINS) { // pt integration
00337             v2IntPt_4[species][cent+1] += v2_4[species][cent+1][ptcount[species]];
00338             usedInt_cent_yields4[species][cent+1] += cent_yields[species][cent][ptbin];
00339           }
00340           v2_4_e[species][cent+1][ptcount[species]] += pow(cent_yields[species][cent][ptbin]*Cumulant_4_e[0][species][cent][ptbin]/(pow(-Cumulant_4[1][species][cent][ptbin],.75)),2);
00341           v2_4_e[species][cent+1][ptcount[species]] += pow(.75*cent_yields[species][cent][ptbin]*Cumulant_4[0][species][cent][ptbin]*Cumulant_4_e[1][species][cent][ptbin]/(pow(-Cumulant_4[1][species][cent][ptbin],1.75)),2);
00342           v2_4_e[species][0][ptcount[species]] += v2_4_e[species][cent+1][ptcount[species]];
00343           if(ptbin<PTINTBINS) v2IntPt_4_e[species][cent+1] += v2_4_e[species][cent+1][ptcount[species]];
00344 
00345           used_pt_yields4[cent+1] = cent_yields[species][cent][ptbin];
00346           used_cent_yields4[species][cent+1] += cent_yields[species][cent][ptbin];
00347           used_pt_yields4[0] += used_pt_yields4[cent+1]; // min bias
00348         }//if
00349         
00350         // 2-particle cumulant v2 
00351         if(Cumulant_2[1][species][cent][ptbin] > 0) { // must be positive
00352           v2_2[species][cent+1][ptcount[species]] += cent_yields[species][cent][ptbin]*Cumulant_2[0][species][cent][ptbin]/(pow(Cumulant_2[1][species][cent][ptbin],.5));
00353           v2_2[species][0][ptcount[species]] += v2_2[species][cent+1][ptcount[species]];
00354           if(ptbin<PTINTBINS) { // pt integration
00355             v2IntPt_2[species][cent+1] += v2_2[species][cent+1][ptcount[species]];
00356             usedInt_cent_yields2[species][cent+1] += cent_yields[species][cent][ptbin];
00357           }       
00358           v2_2_e[species][cent+1][ptcount[species]] += pow(cent_yields[species][cent][ptbin]*Cumulant_2_e[0][species][cent][ptbin]/(pow(Cumulant_2[1][species][cent][ptbin],.5)),2);
00359           v2_2_e[species][cent+1][ptcount[species]] += pow(.5*cent_yields[species][cent][ptbin]*Cumulant_2[0][species][cent][ptbin]*Cumulant_2_e[1][species][cent][ptbin]/(pow(Cumulant_2[1][species][cent][ptbin],1.5)),2);
00360           v2_2_e[species][0][ptcount[species]] += v2_2_e[species][cent+1][ptcount[species]];
00361           if(ptbin<PTINTBINS) v2IntPt_2_e[species][cent+1] += v2_2_e[species][cent+1][ptcount[species]];
00362           
00363           used_pt_yields2[cent+1] = cent_yields[species][cent][ptbin];
00364           used_cent_yields2[species][cent+1] += cent_yields[species][cent][ptbin];
00365           used_pt_yields2[0] += used_pt_yields2[cent+1]; // min bias
00366         }//if
00367                                     
00368       }//cent
00369       
00370       for(int cent=0; cent<CENTBINS+1; cent++){ // includes cent=0
00371         
00372         if(used_pt_yields2[cent] > 0){
00373           v2_2[species][cent][ptcount[species]] /= used_pt_yields2[cent];
00374           v2_2_e[species][cent][ptcount[species]] = sqrt(v2_2_e[species][cent][ptcount[species]]);
00375           v2_2_e[species][cent][ptcount[species]] /= used_pt_yields2[cent];
00376         }
00377         if(used_pt_yields4[cent] > 0){
00378           v2_4[species][cent][ptcount[species]] /= used_pt_yields4[cent];
00379           v2_4_e[species][cent][ptcount[species]] = sqrt(v2_4_e[species][cent][ptcount[species]]);
00380           v2_4_e[species][cent][ptcount[species]] /= used_pt_yields4[cent];
00381         }
00382 
00383       }//cent
00384       
00385       ptcount[species]++;
00386     }//pt
00387 
00388     // For pt integrated v2
00389     for(int cent=0; cent<CENTBINS+1; cent++){ // includes cent=0
00390       if(usedInt_cent_yields2[species][cent] > 0){
00391         v2IntPt_2[species][cent] /= usedInt_cent_yields2[species][cent];
00392         v2IntPt_2_e[species][cent] = sqrt(v2IntPt_2_e[species][cent]);
00393         v2IntPt_2_e[species][cent] /= usedInt_cent_yields2[species][cent];
00394       }
00395       if(usedInt_cent_yields4[species][cent] > 0){
00396         v2IntPt_4[species][cent] /= usedInt_cent_yields4[species][cent];
00397         v2IntPt_4_e[species][cent] = sqrt(v2IntPt_4_e[species][cent]);
00398         v2IntPt_4_e[species][cent] /= usedInt_cent_yields4[species][cent];
00399       }
00400     }//cent
00401 
00402   }//species
00403   
00404   cout << endl;
00405   for(int ptbin=0; ptbin<PTBINS; ptbin++){ // for checking
00406     //cout << ptbin << ": " << v2_2[0][3][ptbin] << " +/- " << v2_2_e[0][3][ptbin] << endl;
00407     //cout << ptbin << ": " << cent_yields[0][3][ptbin] << endl;
00408     //cout << ptbin << ": " << Cumulant_2[1][0][5][ptbin] << endl;
00409   }
00410 
00411   cout << endl << "cent: \t v2IntPt_2 +/- err,\t v2IntPt_4 +/- err" << endl;
00412   for(int cent=0; cent<CENTBINS+1; cent++){
00413 
00414     //Throw away last point since its middle-pt value is not well defined.
00415     //It represents all particles with pt > 6 GeV/c
00416     //cout << ptOverflow << ": " << v2_2[0][cent][ptOverflow] << endl;
00417     v2_2[0][cent][ptOverflow]= 10.;
00418     v2_4[0][cent][ptOverflow]= 10.;
00419     //Throw away next-to-last point since it contains the integrated cumulant.
00420     v2_2[0][cent][ptCu]= 10.;
00421     v2_4[0][cent][ptCu]= 10.;
00422     //cout << ptOverflow << ": " << v2_2[0][cent][ptOverflow] << endl; // pt = 10?
00423 
00424     // output integrated over pt
00425     cout << cent << ": \t" << setprecision(3) << v2IntPt_2[0][cent]*100. << " \t+/- " <<
00426       v2IntPt_2_e[0][cent]*100. << ", \t" << v2IntPt_4[0][cent]*100. << " \t+/- " <<
00427       v2IntPt_4_e[0][cent]*100. << endl;
00428 
00429   }//cent
00430   cout << endl;
00431 
00433   // plot the results 
00434 
00435   TFile graphFile("flow.dirCumulant.graphs.root", "RECREATE"); // output
00436 
00437   // For v2(pt)
00438 
00439   double middle_pt_points[PTBINS]={0};
00440   for(int i=0; i<PTBINS; i++){
00441     middle_pt_points[i]=double((i)/10. + 0.1); // GeV = bin/10.
00442   }
00443   double middle_pt_points_e[PTBINS]={0}; // no errors
00444   
00445   float v2max = 0.4;
00446   
00447   TGraphErrors *charged_v2_2[CENTBINS+1];
00448   TGraphErrors *charged_v2_4[CENTBINS+1];
00449   TCanvas* can[CENTBINS+1];
00450   int canvasWidth = 780, canvasHeight = 600;             // landscape
00451   for(int cent=0; cent<CENTBINS+1; cent++){
00452     TString* canName = new TString("centrality_"); 
00453     *canName += (cent); 
00454     if(ptGraphs) {
00455       cout << "plot " << canName->Data() << endl;
00456       can[cent] = new TCanvas(canName->Data(), canName->Data(), canvasWidth, canvasHeight);
00457     }
00458     TLegend *legend = new TLegend(.1,.7,.46,.9,NULL,"brNDC");
00459     legend->SetFillColor(kWhite);
00460     TPaveLabel* title = new TPaveLabel(0.7,0.96,0.9,0.99,canName->Data());
00461     if(ptGraphs) title->Draw();
00462 
00463     TString graphName("v22_"); 
00464     graphName += (cent); 
00465     charged_v2_2[cent]= new TGraphErrors(PTBINS,middle_pt_points,v2_2[0][cent],middle_pt_points_e,v2_2_e[0][cent]);
00466     charged_v2_2[cent]->SetMarkerStyle(20);
00467     charged_v2_2[cent]->SetMarkerColor(2);
00468     charged_v2_2[cent]->SetLineColor(2);
00469     charged_v2_2[cent]->SetMinimum(0.);
00470     charged_v2_2[cent]->SetMaximum(v2max);
00471     charged_v2_2[cent]->GetXaxis()->SetTitle("P_{t} (GeV/c)");
00472     charged_v2_2[cent]->GetYaxis()->SetTitle("v_{2}");
00473     charged_v2_2[cent]->SetTitle("Elliptic Flow");
00474     charged_v2_2[cent]->Write(graphName.Data());
00475     
00476     TString graphName("v24_"); 
00477     graphName += (cent); 
00478     charged_v2_4[cent]= new TGraphErrors(PTBINS,middle_pt_points,v2_4[0][cent],middle_pt_points_e,v2_4_e[0][cent]);
00479     charged_v2_4[cent]->SetMarkerStyle(20);
00480     charged_v2_4[cent]->SetMarkerColor(4);
00481     charged_v2_4[cent]->SetLineColor(4);
00482     charged_v2_4[cent]->SetMinimum(0.);
00483     charged_v2_4[cent]->SetMaximum(v2max);
00484     charged_v2_4[cent]->GetXaxis()->SetTitle("P_{t} (GeV/c)");
00485     charged_v2_4[cent]->GetYaxis()->SetTitle("v_{2}{4}");
00486     charged_v2_4[cent]->SetTitle("Elliptic Flow");
00487     charged_v2_4[cent]->Write(graphName.Data());
00488     
00489     if(ptGraphs) charged_v2_2[cent]->Draw("AP");
00490     legend->AddEntry(charged_v2_2[cent],"charged hadron v_{2}{2}","p");
00491     if(ptGraphs) charged_v2_4[cent]->Draw("P");
00492     legend->AddEntry(charged_v2_4[cent],"charged hadron v_{2}{4}","p");    
00493     if(ptGraphs) legend->Draw("same");
00494     
00495     //can[cent]->Print(".pdf");
00496 
00497   }
00498 
00499   //Event_counter histograms
00500 
00501   TH1* Event_counter = dynamic_cast<TH1*>(myfile->Get("Event_counter"));
00502   if (!Event_counter) {
00503     cout << "### Can't find hist Event_counter" << endl;
00504   } else {
00505     // make the graph page
00506     gStyle->SetOptStat("e");
00507     TString* canEvtName = new TString("Event_counter"); 
00508     cout << "plot " << canEvtName->Data() << endl;
00509     TCanvas* canEvt = new TCanvas(canEvtName->Data(), canEvtName->Data(), 780, 600);
00510     Event_counter->SetMinimum(0.);
00511     Event_counter->GetXaxis()->SetTitle("centrality bin");
00512     Event_counter->Draw();
00513     Event_counter->Write(canEvtName->Data());    
00514   }
00515   TH1* Event_counterWeighted = dynamic_cast<TH1*>(myfile->Get("Event_counterWeighted"));
00516   if (!Event_counterWeighted) {
00517     cout << "### Can't find hist Event_counterWeighted" << endl;
00518   } else {
00519     // make the graph page
00520     //gStyle->SetOptStat("e");
00521     TString* canEvtWgtName = new TString("Event_counterWeighted"); 
00522     cout << "plot " << canEvtWgtName->Data() << endl;
00523     TCanvas* canEvtWgt = new TCanvas(canEvtWgtName->Data(), canEvtWgtName->Data(), 780, 600);
00524     Event_counterWeighted->SetMinimum(0.);
00525     Event_counterWeighted->GetXaxis()->SetTitle("centrality bin");
00526     Event_counterWeighted->Draw("hist");
00527     Event_counterWeighted->Write(canEvtWgtName->Data());    
00528   }
00529   //gStyle->SetOptStat(0);
00530 
00531   // v2 integrated
00532   
00533   float v2Intmax = 8.;
00534 
00535   // make the graph page
00536   TString* canIntName = new TString("v2Int"); 
00537   cout << "plot " << canIntName->Data() << endl;
00538   TCanvas* canInt = new TCanvas(canIntName->Data(), canIntName->Data(), 600, 780);
00539   TLegend* legendInt = new TLegend(0.35,0.1,0.6,0.35);
00540   legendInt->SetFillColor(kWhite);
00541 
00542   // make a histogram
00543   TString* histCenName = new TString("v2Int");
00544   TH1F* histCen = new TH1F(histCenName->Data(), histCenName->Data(), 80, 0., 80.);
00545   histCen->SetLineColor(kBlack);
00546   histCen->Draw();
00547 
00548   TGraphErrors *grv2IntPt_2 = new TGraphErrors(CENTBINS);
00549   TGraphErrors *grv2IntPt_4 = new TGraphErrors(CENTBINS);
00550   TGraphErrors *grv2IntCu_2 = new TGraphErrors(CENTBINS);
00551   TGraphErrors *grv2IntCu_4 = new TGraphErrors(CENTBINS);
00552 
00553   int m = 0;
00554   for(int cent=0; cent<CENTBINS; cent++){
00555     grv2IntPt_2->SetPoint(m, x[CENTBINS-cent-1], v2IntPt_2[0][cent+1]*100.);
00556     grv2IntPt_2->SetPointError(m, 0., v2IntPt_2_e[0][cent+1]*100.);
00557     grv2IntPt_4->SetPoint(m, x[CENTBINS-cent-1], v2IntPt_4[0][cent+1]*100.);
00558     grv2IntPt_4->SetPointError(m, 0., v2IntPt_4_e[0][cent+1]*100.);
00559     grv2IntCu_2->SetPoint(m, x[CENTBINS-cent-1], v2IntCu_2[cent]*100.);
00560     grv2IntCu_2->SetPointError(m, 0., v2IntCu_2Err[cent]*100.);
00561     grv2IntCu_4->SetPoint(m, x[CENTBINS-cent-1], v2IntCu_4[cent]*100.);
00562     grv2IntCu_4->SetPointError(m, 0., v2IntCu_4Err[cent]*100.);
00563 
00564     m++;
00565   }//cent
00566 
00567   grv2IntPt_2->SetMarkerStyle(kFullCircle);
00568   grv2IntPt_2->SetMarkerColor(kRed);
00569   grv2IntPt_2->SetLineColor(kRed);
00570   grv2IntPt_2->SetMinimum(0.);
00571   grv2IntPt_2->SetMaximum(v2Intmax);
00572   grv2IntPt_2->SetTitle("Elliptic Flow");
00573   grv2IntPt_2->Draw("AP");
00574   legendInt->AddEntry(grv2IntPt_2,"v_{2}{2}(pt)","p");
00575 
00576   grv2IntPt_4->SetMarkerStyle(kFullSquare);
00577   grv2IntPt_4->SetMarkerColor(kBlue);
00578   grv2IntPt_4->SetLineColor(kBlue);
00579   grv2IntPt_4->Draw("P");
00580   legendInt->AddEntry(grv2IntPt_4,"v_{2}{4}(pt)","p");
00581 
00582   grv2IntCu_2->SetMarkerStyle(kOpenCircle);
00583   grv2IntCu_2->SetMarkerColor(kRed);
00584   grv2IntCu_2->SetLineColor(kRed);
00585   grv2IntCu_2->Draw("P");
00586   legendInt->AddEntry(grv2IntCu_2,"v_{2}{2}(cu)","p");
00587 
00588   grv2IntCu_4->SetMarkerStyle(kOpenSquare);
00589   grv2IntCu_4->SetMarkerColor(kBlue);
00590   grv2IntCu_4->SetLineColor(kBlue);
00591   grv2IntCu_4->Draw("P");
00592   legendInt->AddEntry(grv2IntCu_4,"v_{2}{4}(cu)","p");
00593     
00594   legendInt->Draw("same");
00595   // output
00596   grv2IntPt_2->Write("v22IntPt");
00597   grv2IntPt_4->Write("v24IntPt");
00598   grv2IntCu_2->Write("v22IntCu");
00599   grv2IntCu_4->Write("v24IntCu");
00600     
00601   //can[cent]->Print(".pdf");
00602       
00603   // label the axes
00604   TLatex l;
00605   l.SetNDC();
00606   l.SetTextSize(0.06);
00607   l.DrawLatex(0.5,0.02,"% Most Central");
00608   l.SetTextAngle(90);
00609   l.DrawLatex(0.05,0.7,"v_{2} (%)" );
00610 
00611   graphFile.Close();
00612   
00613 }
StFlowAnalysisMaker/directCumulants_v2.C
