StRoot: StFlowAnalysisMaker/calculateCumulant.C Source File

00001 #include "float.h"
00002 #include "calculateCumulant.h"
00003 
00004 void calculateCumulant(const char* histFileName){
00005 
00006   //  gSystem->Load("myProjectProfile_C.so"); //user defined integration
00007 
00008   double r0Sq=r0*r0;
00009   double perCent=0.01;
00010 
00011 
00012   TString xLabel = "Pseudorapidity";
00013   if (isPidFlow) { xLabel = "Rapidity"; }  
00014   
00015 
00016   TFile* f = new TFile(histFileName,"UPDATE");
00017 
00018   struct histFullHars {
00019 
00020     TProfile2D**   mHistCumul2D;
00021     TProfile**     mHistCumulEta;
00022     TProfile**     mHistCumulPt;
00023 
00024     TProfile2D*    mHistCumulMix2D;
00025     TProfile*      mHistCumulMixEta;
00026     TProfile*      mHistCumulMixPt;
00027 
00028     TH2D**         mHist_v2D;
00029     TH1D**         mHist_vEta;
00030     TH1D**         mHist_vPt;
00031 
00032 
00033 
00034     TH2D*         mHistMix_v2D;
00035     TH1D*         mHistMix_vEta;
00036     TH1D*         mHistMix_vPt;
00037 
00038 
00039 
00040 
00041     Double_t       mCumulIntegG0[nCumulIntegOrders*nCumulInteg_qMax]; 
00042     TH1D*          mHistMultSum;    //histo for holding mMultSum.
00043     TH1D*          mHistNEvent;
00044 
00045     TH1D*          mHistMeanWgtSqrSum; // sum <w^2> over evts.
00046 
00047     Double_t       mCumulIntegG0Mix[nCumulMixHar_pMax*nCumulMixHar_qMax];  
00048 
00049     TH1D**         mHistCumulIntegG0Sum; //summation of G value.
00050     TH1D**         mHistCumulIntegG0MixSum; //summation of G value.
00051 
00052   };
00053 
00054   struct histFulls {
00055    TProfile**  mHistCumul; //integrated cumulant from differential
00056    TH1D**      mHist_v; //integrated flow from differential
00057 
00058    TProfile*   mHistCumulMix; //cumulant for mix 1st and 2nd hars.
00059    TH1D*       mHistMix_v; //integrated flow from differential
00060 
00061     struct histFullHars histFullHar[nHars];
00062   };
00063 
00064   struct histFulls histFull[nSels]; 
00065 
00066 
00067   TString* histTitle;
00068   TString* histTitle2;
00069 
00070 
00071   //================        fill in pointers    ===================
00072   //for each selection
00073   for (int k = 0; k < nSels; k++) {
00074     char countSels[2];
00075     sprintf(countSels,"%d",k+1);
00076 
00077     histFull[k].mHistCumul = new TProfile*[nCumulDiffOrders];
00078     histFull[k].mHist_v    = new TH1D*[nCumulDiffOrders];
00079   
00080     for (int ord = 0; ord < nCumulDiffOrders; ord++) {
00081       char theCumulOrderChar[2]; // if >10, need to use char*
00082       sprintf(theCumulOrderChar,"%d",(ord+1)*2);
00083       histTitle = new TString("Flow_Cumul_Order"); 
00084       histTitle->Append(*theCumulOrderChar);           
00085       histTitle->Append("_Sel");                      
00086       histTitle->Append(*countSels);        
00087       histFull[k].mHistCumul[ord] = (TProfile *)f->Get(histTitle->Data());
00088       delete histTitle;
00089 
00090       histTitle = new TString("Flow_Cumul_v_Order");
00091       histTitle->Append(*theCumulOrderChar);
00092       histTitle->Append("_Sel");
00093       histTitle->Append(*countSels);
00094       histFull[k].mHist_v[ord] = 
00095         new  TH1D(*(histFull[k].mHistCumul[ord]->ProjectionX(histTitle->Data(),"e")));
00096       histFull[k].mHist_v[ord]->SetTitle(histTitle->Data());
00097       histFull[k].mHist_v[ord]->SetXTitle("harmonic");
00098       histFull[k].mHist_v[ord]->SetYTitle("v (%)");
00099       if (ord>0) histFull[k].mHist_v[ord]->Scale(1./profScale);
00100       delete histTitle;
00101 
00102     }
00103 
00104 
00106     
00107       histTitle = new TString("Flow_CumulMix");
00108       histTitle->Append("_Sel");                      
00109       histTitle->Append(*countSels); 
00110       histFull[k].mHistCumulMix = (TProfile *)f->Get(histTitle->Data());
00111       delete histTitle;
00112 
00113 
00114       histTitle = new TString("Flow_CumulMix_v_Sel");
00115       histTitle->Append(*countSels);
00116       histFull[k].mHistMix_v = 
00117         new  TH1D(*(histFull[k].mHistCumulMix->ProjectionX(histTitle->Data(),"e")));
00118       histFull[k].mHistMix_v->SetTitle(histTitle->Data());
00119       histFull[k].mHistMix_v->SetXTitle("place for v1 from mixed harmonic");
00120       histFull[k].mHistMix_v->SetYTitle("v (%)");
00121       histFull[k].mHistMix_v->Scale(1./profScale);
00122       delete histTitle;
00124 
00125 
00126 
00127 
00128     // for each harmonic
00129     for (int j = 0; j < nHars; j++) {
00130       char countHars[2];
00131       sprintf(countHars,"%d",j+1);
00132 
00133       histFull[k].histFullHar[j].mHistCumul2D  =  
00134         new TProfile2D*[nCumulDiffOrders];
00135       histFull[k].histFullHar[j].mHistCumulEta = 
00136         new TProfile*[nCumulDiffOrders];
00137       histFull[k].histFullHar[j].mHistCumulPt  = 
00138         new TProfile*[nCumulDiffOrders];
00139 
00140       histFull[k].histFullHar[j].mHist_v2D  = 
00141         new TH2D*[nCumulDiffOrders];
00142       histFull[k].histFullHar[j].mHist_vEta = 
00143         new TH1D*[nCumulDiffOrders];
00144       histFull[k].histFullHar[j].mHist_vPt  = 
00145         new TH1D*[nCumulDiffOrders];
00146 
00147 
00149       if (j==0){ //only j==0 makes sense here.
00150 
00151         histTitle = new TString("Flow_CumulMix2D");
00152         histTitle->Append("_Sel");
00153         histTitle->Append(*countSels);
00154         histTitle->Append("_Har");
00155         histTitle->Append(*countHars);
00156         histFull[k].histFullHar[j].mHistCumulMix2D  = 
00157                (TProfile2D *)f->Get(histTitle->Data());
00158         delete histTitle;
00159 
00160 
00161 
00162         histTitle = new TString("Flow_CumulMixEta");
00163         histTitle->Append("_Sel");
00164         histTitle->Append(*countSels);
00165         histTitle->Append("_Har");
00166         histTitle->Append(*countHars);
00167         histFull[k].histFullHar[j].mHistCumulMixEta =  
00168                (TProfile *)f->Get(histTitle->Data());
00169         delete histTitle;
00170 
00171         histTitle = new TString("Flow_CumulMixPt");
00172         histTitle->Append("_Sel");
00173         histTitle->Append(*countSels);
00174         histTitle->Append("_Har");
00175         histTitle->Append(*countHars);
00176         histFull[k].histFullHar[j].mHistCumulMixPt =  
00177                (TProfile *)f->Get(histTitle->Data());
00178         delete histTitle;
00179 
00180 
00181 
00182       histFull[k].histFullHar[j].mHistCumulIntegG0MixSum =
00183         new       TH1D*[nCumulMixHar_pMax*nCumulMixHar_qMax];
00184 
00185 
00186     for (int pq  = 0; pq < nCumulMixHar_pMax*nCumulMixHar_qMax ; pq++) { 
00187         int pIndex = (pq/nCumulMixHar_qMax);
00188         int qIndex        = pq%nCumulMixHar_qMax;
00189         
00190         char pIndexChar[2];
00191         char qIndexChar[2]; // if >10, need to use char*
00192         sprintf(pIndexChar,"%d",pIndex); 
00193         sprintf(qIndexChar,"%d",qIndex);
00194 
00195         histTitle = new TString("Flow_CumulIntegG0MixSum");
00196         histTitle->Append("_GenFunIdxp");
00197         histTitle->Append(*pIndexChar);
00198         histTitle->Append("_GenFunIdxq");
00199         histTitle->Append(*qIndexChar);
00200         histTitle->Append("_Sel");
00201         histTitle->Append(*countSels);
00202         histTitle->Append("_Har");
00203         histTitle->Append(*countHars);
00204         histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq]=
00205                (TH1D *)f->Get(histTitle->Data());
00206         delete histTitle;  
00207 
00208         histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] = 0.;
00209     }
00210       }
00211 
00212 
00214 
00215 
00216 
00217       // For each cumulant order
00218       for (int ord = 0; ord < nCumulDiffOrders; ord++) {
00219         char theCumulOrderChar[2]; // if >10, need to use char*
00220         sprintf(theCumulOrderChar,"%d",(ord+1)*2);
00221         
00222         histTitle = new TString("Flow_Cumul2D_Order");
00223         histTitle->Append(*theCumulOrderChar);
00224         histTitle->Append("_Sel");
00225         histTitle->Append(*countSels);
00226         histTitle->Append("_Har");
00227         histTitle->Append(*countHars);
00228         histFull[k].histFullHar[j].mHistCumul2D[ord] = 
00229                (TProfile2D *)f->Get(histTitle->Data());
00230         delete histTitle;
00231 
00232         histTitle = new TString("Flow_CumulEta_Order");
00233         histTitle->Append(*theCumulOrderChar);
00234         histTitle->Append("_Sel");
00235         histTitle->Append(*countSels);
00236         histTitle->Append("_Har");
00237         histTitle->Append(*countHars);
00238         histFull[k].histFullHar[j].mHistCumulEta[ord]= 
00239                (TProfile *)f->Get(histTitle->Data());
00240         delete histTitle;
00241 
00242 
00243         histTitle = new TString("Flow_CumulPt_Order");
00244         histTitle->Append(*theCumulOrderChar);
00245         histTitle->Append("_Sel");
00246         histTitle->Append(*countSels);
00247         histTitle->Append("_Har");
00248         histTitle->Append(*countHars);
00249         histFull[k].histFullHar[j].mHistCumulPt[ord] = 
00250                (TProfile *)f->Get(histTitle->Data());
00251         delete histTitle;  
00252 
00253       }
00254 
00255 
00256       histFull[k].histFullHar[j].mHistCumulIntegG0Sum =
00257       new TH1D*[nCumulIntegOrders*nCumulInteg_qMax];
00258 
00259       for (int pq = 0; pq <  nCumulIntegOrders*nCumulInteg_qMax; pq++) {
00260         int cumulIndex = (pq/nCumulInteg_qMax) + 1; // like 1,2,3. 
00261         //not begining with 0. That is "p" in Eq. (B1, PG5).
00262         int qIndex = pq%(nCumulInteg_qMax); // like 0,1,..qMax-1 
00263         //begining with 0. Just like Eq. (B3, PG5).
00264 
00265         char theCumulOrderChar[2];
00266         char qIndexOrderChar[2]; // if >10, need to use char*
00267         sprintf(theCumulOrderChar,"%d",cumulIndex*2); 
00268         sprintf(qIndexOrderChar,"%d",qIndex);
00269 
00270         histTitle = new TString("Flow_CumulIntegG0Sum_Order");
00271         histTitle->Append(*theCumulOrderChar);
00272         histTitle->Append("_GenFunIdx");
00273         histTitle->Append(*qIndexOrderChar);
00274         histTitle->Append("_Sel");
00275         histTitle->Append(*countSels);
00276         histTitle->Append("_Har");
00277         histTitle->Append(*countHars);
00278         histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq] =
00279                (TH1D *)f->Get(histTitle->Data());
00280         delete histTitle;  
00281 
00282         histFull[k].histFullHar[j].mCumulIntegG0[pq] = 0.;
00283 
00284       }
00285 
00286 
00287 
00288 
00289       histTitle = new TString("Flow_CumulMultSum_Sel");
00290       histTitle->Append(*countSels);
00291       histTitle->Append("_Har");
00292       histTitle->Append(*countHars);
00293       histFull[k].histFullHar[j].mHistMultSum =
00294              (TH1D *)f->Get(histTitle->Data());
00295         delete histTitle;  
00296 
00297 
00298       histTitle = new TString("Flow_CumulNEvent_Sel");
00299       histTitle->Append(*countSels);
00300       histTitle->Append("_Har");
00301       histTitle->Append(*countHars);
00302       histFull[k].histFullHar[j].mHistNEvent =
00303              (TH1D *)f->Get(histTitle->Data());
00304         delete histTitle;  
00305 
00306 
00307       histTitle = new TString("Flow_CumulMeanWgtSqrSum_Sel");
00308       histTitle->Append(*countSels);
00309       histTitle->Append("_Har");
00310       histTitle->Append(*countHars);
00311       histFull[k].histFullHar[j].mHistMeanWgtSqrSum =
00312              (TH1D *)f->Get(histTitle->Data());
00313       delete histTitle;
00314 
00315 
00316     }
00317   }
00318 
00319  //================   end of fill in pointers    ===================
00320 
00321 
00322 
00323 
00324   for (int k = 0; k < nSels; k++) {
00325 
00326     char countSels[2];
00327     sprintf(countSels,"%d",k+1);
00328 
00329     double  meanIntegV[nHars];     // V**1
00330     double  meanIntegV2[nHars];    // V**2
00331     double  meanIntegV3[nHars];    // V**3
00332     double  meanIntegV4[nHars];    // V**4
00333     double  cumulInteg1[nHars];    // outside of harmonic loop
00334     double  cumulInteg2[nHars];
00335     double  cumulInteg3[nHars];
00336     double  q2[nHars]; // for old method. <Q>**2 in (74) of old paper.
00337     double  q4[nHars];
00338     double  q6[nHars];
00339 
00340 
00341     double cumulIntegMix[nHars];
00342     double meanIntegVMix[nHars];
00343     
00344     for (int j = 0; j < nHars; j++) {
00345       meanIntegV[j]  = 0.;
00346       meanIntegV2[j] = 0.;
00347       meanIntegV3[j] = 0.;
00348       meanIntegV4[j] = 0.;
00349       cumulInteg1[j] = 0.;
00350       cumulInteg2[j] = 0.;
00351       cumulInteg3[j] = 0.;
00352       cumulIntegMix[j] = 0.;
00353       meanIntegVMix[j] = 0.;
00354     }
00355 
00356     for (int j = 0; j < nHars; j++) {
00357 
00358 
00359       char countHars[2];
00360       sprintf(countHars,"%d",j+1);
00361 
00362       cout<<" ========== Sel"<<k+1<<" Har"<<j+1<<" ==========="<<endl;
00363 
00364       cout<<" event # "<<histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1)<<endl;
00365 
00366       double mAvMult = 
00367         histFull[k].histFullHar[j].mHistMultSum->GetBinContent(1)/
00368         histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
00369 
00370       cout<<"mAvMult Sel"<<k<<" har"<<j<<" "<<mAvMult<<endl;
00371 
00372 
00373       double CpInteg[nCumulIntegOrders]; // Cp in (B4, PG6)
00374       
00375       for (int pq = 0; pq < nCumulIntegOrders; pq ++) CpInteg[pq] = 0.;
00376       for (int pq = 0; pq < nCumulIntegOrders*nCumulInteg_qMax; pq++) {   
00377 
00378         int theCumulOrder = (pq/nCumulInteg_qMax) + 1; // like 1,2,3.   
00379 
00380         histFull[k].histFullHar[j].mCumulIntegG0[pq] =
00381         histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq]->GetBinContent(1)/
00382         histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
00383 
00384         //      cout<<" --- histFull[k].histFullHar[j].mCumulIntegG0[pq] "<<histFull[k].histFullHar[j].mCumulIntegG0[pq]<<endl;
00385 
00386         if (!isNewMethod){
00387           CpInteg[theCumulOrder-1] +=
00388             (log(histFull[k].histFullHar[j].mCumulIntegG0[pq]) /
00389              ((float)nCumulInteg_qMax));
00390         } else {
00391           CpInteg[theCumulOrder-1] +=
00392             (mAvMult*(pow(histFull[k].histFullHar[j].mCumulIntegG0[pq], 1./mAvMult) -1.) /
00393              float(nCumulInteg_qMax)); // (B3, PG6) 
00394           //          cout<<" CpInteg[theCumulOrder-1] "<<CpInteg[theCumulOrder-1]<<endl;
00395         }
00396       }
00397 
00398       cumulInteg1[j] = // (B5, PG7)
00399         (3.*CpInteg[1-1] - (3./2.)*CpInteg[2-1] + (1./3.)*CpInteg[3-1]) / r0Sq;
00400       
00401       cumulInteg2[j] = ((-10.*CpInteg[1-1]) + (8.*CpInteg[2-1]) - 
00402                          (2.*CpInteg[3-1])) / (r0Sq*r0Sq); 
00403       
00404       cumulInteg3[j] = ( (18.*CpInteg[1-1]) - (18.*CpInteg[2-1]) + (6.*CpInteg[3-1]))
00405                         / (r0Sq*r0Sq*r0Sq);
00406       
00407 
00408       for (int ord = 0; ord < nCumulDiffOrders; ord++) {
00409         char theCumulOrderChar[2]; // if >10, need to use char*
00410         sprintf(theCumulOrderChar,"%d",(ord+1)*2);
00411         
00412 
00413         
00414         histTitle = new TString("Flow_Cumul_vEta_Order");
00415         histTitle->Append(*theCumulOrderChar);
00416         histTitle->Append("_Sel");
00417         histTitle->Append(*countSels);
00418         histTitle->Append("_Har");
00419         histTitle->Append(*countHars);
00420         histFull[k].histFullHar[j].mHist_vEta[ord] = 
00421           new  TH1D(*(histFull[k].histFullHar[j].mHistCumulEta[ord]->
00422                       ProjectionX(histTitle->Data(),"e")));
00423         histFull[k].histFullHar[j].mHist_vEta[ord]->SetTitle(histTitle->Data());
00424         histFull[k].histFullHar[j].mHist_vEta[ord]->SetXTitle((char*)xLabel.Data());
00425         histFull[k].histFullHar[j].mHist_vEta[ord]->SetYTitle("v (%)");
00426         if (ord>0) histFull[k].histFullHar[j].mHist_vEta[ord]->Scale(1./profScale);
00427         delete histTitle;
00428         
00429         histTitle = new TString("Flow_Cumul_vPt_Order");
00430         histTitle->Append(*theCumulOrderChar);
00431         histTitle->Append("_Sel");
00432         histTitle->Append(*countSels);
00433         histTitle->Append("_Har");
00434         histTitle->Append(*countHars);
00435 
00436         //    if you want integrate for only TPC/FTPC, you need to do
00437         //  1. load myProjectProfile.C at the begining
00438         //  2. append "_TPC" or "_FTPC" here to the title
00439         //  3. use myProjectionY (X) to get mHist_vPt or mHist_vEta.
00440         //              histTitle->Append("_TPC");
00441 
00442           histTitle2 = new TString("Flow_Cumul2D_Order");
00443           histTitle2->Append(*theCumulOrderChar);
00444           histTitle2->Append("_Sel");
00445           histTitle2->Append(*countSels);
00446           histTitle2->Append("_Har");
00447           histTitle2->Append(*countHars);
00448 
00449           histFull[k].histFullHar[j].mHist_vPt[ord] = 
00450             //uncomment this line if you want to integrate over different range
00451             //myProjectionY((TProfile2D *)f->Get(histTitle2->Data()) , 0., 1., 1, ((j==0 || j== 2) ? 1 : 0) );
00452                           new  TH1D(*(histFull[k].histFullHar[j].mHistCumulPt[ord]->ProjectionX(histTitle->Data(),"e")));
00453         histFull[k].histFullHar[j].mHist_vPt[ord]->SetTitle(histTitle->Data());
00454         histFull[k].histFullHar[j].mHist_vPt[ord]->SetName(histTitle->Data());
00455         histFull[k].histFullHar[j].mHist_vPt[ord]->SetXTitle("Pt (GeV)");
00456         histFull[k].histFullHar[j].mHist_vPt[ord]->SetYTitle("v (%)");
00457         if (ord>0) histFull[k].histFullHar[j].mHist_vPt[ord]->Scale(1./profScale);
00458         delete histTitle; 
00459         delete histTitle2;
00460 
00461 
00462         //do not reverse the order of getting v2D and vPt !!
00463 
00464         histTitle = new TString("Flow_Cumul_v2D_Order");
00465         histTitle->Append(*theCumulOrderChar);
00466         histTitle->Append("_Sel");
00467         histTitle->Append(*countSels);
00468         histTitle->Append("_Har");
00469         histTitle->Append(*countHars);
00470         histFull[k].histFullHar[j].mHist_v2D[ord] = 
00471           new TH2D(*(histFull[k].histFullHar[j].mHistCumul2D[ord]->
00472                      ProjectionXY(histTitle->Data(),"e")));
00473         histFull[k].histFullHar[j].mHist_v2D[ord]->SetTitle(histTitle->Data());
00474         histFull[k].histFullHar[j].mHist_v2D[ord]->SetXTitle((char*)xLabel.Data());
00475         histFull[k].histFullHar[j].mHist_v2D[ord]->SetYTitle("Pt (GeV)");
00476         histFull[k].histFullHar[j].mHist_v2D[ord]->SetZTitle("v (%)");
00477         if (ord>0) histFull[k].histFullHar[j].mHist_v2D[ord]->Scale(1./profScale);
00478         delete histTitle;
00479 
00480 
00481 
00482       }
00483 
00484 
00485       //In the case if cumulant fails, let <V> etc. to be 1. to avoid crash.
00486         meanIntegV[j]  = (cumulInteg1[j]>0.) ? sqrt(cumulInteg1[j]) : 1.;           // <v>    2-part, m=1
00487         meanIntegV2[j] = (cumulInteg1[j]>0.) ? cumulInteg1[j] : 1.;                 // <v**2> 2-part, m=2
00488         meanIntegV3[j] = (cumulInteg2[j]<0.) ? pow(-1.*cumulInteg2[j], 3./4.) : 1.; // <v**3> 4-part, m=1
00489         meanIntegV4[j] = (cumulInteg2[j]<0.) ? -1.*cumulInteg2[j] : 1.;             // <v**4> 4-part, m=2
00490 
00491 
00492     
00493 
00494       if (m_M==1) { // Eq. (PG14) 
00495 
00496         histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
00497         histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-1./(meanIntegV3[j]*perCent));// (34b)
00498         histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
00499         histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
00500         histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
00501         histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
00502       } else if (m_M==2) {
00503         histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
00504         histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
00505         histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV2[j]*perCent));   // (35a)
00506         histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-0.5/(meanIntegV4[j]*perCent) ); // (35b)
00507         histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
00508         histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
00509       }
00510 
00511         if (cumulInteg1[j]<0.) {
00512          cout<<" Sel"<<k+1<<", <V**2> less than zero ! v"<<j+1<<" from 2 particle correlation failed."<<endl;
00513          histFull[k].histFullHar[j].mHist_v2D[0]->Reset();
00514          histFull[k].histFullHar[j].mHist_vEta[0]->Reset();
00515          histFull[k].histFullHar[j].mHist_vPt[0]->Reset();
00516         }
00517 
00518 
00519         if (-1.*cumulInteg2[j]<0.) {
00520          cout<<" Sel"<<k+1<<", <V**4> less than zero ! v"<<j+1<<" from 4 particle correlation failed."<<endl;
00521          histFull[k].histFullHar[j].mHist_v2D[1]->Reset();
00522          histFull[k].histFullHar[j].mHist_vEta[1]->Reset();
00523          histFull[k].histFullHar[j].mHist_vPt[1]->Reset();
00524         }
00525 
00526     }// end of  for (int j = 0; j < nHars; j++) 
00527 
00528 
00529 
00530 
00532 
00533 
00534     for (int j = 0; j < nHars; j++) {
00535 
00536       if (j != 0) continue; //only j==0 makes sense for mixed har.
00537           
00538       char countHars[2];
00539       sprintf(countHars,"%d",j+1);
00540       
00541       cout<<" ========== Sel"<<k+1<<" Har"<<j+1<<" ===== v1{3} ======"<<endl;
00542 
00543       double mAvMult = 
00544         histFull[k].histFullHar[j].mHistMultSum->GetBinContent(1)/
00545         histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
00546 
00547       //      cout<<" mAveMult "<<mAvMult<<endl;
00548 
00549       double mAveMeanWgtSqr = // <w**2>
00550         histFull[k].histFullHar[j].mHistMeanWgtSqrSum->GetBinContent(1)/
00551         histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
00552 
00553       cout<<" histFull[k].histFullHar[j].mHistMeanWgtSqrSum->GetBinContent(1) "<<histFull[k].histFullHar[j].mHistMeanWgtSqrSum->GetBinContent(1)<<endl;
00554       cout<<" histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1) "<<histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1)<<endl;
00555 
00556 
00557        double  CpqMix[nCumulMixHar_pMax][nCumulMixHar_qMax];//(30)
00558           for (int pq = 0; pq < nCumulMixHar_pMax*nCumulMixHar_qMax; pq++) {
00559                 int pIndex = pq/nCumulMixHar_qMax;
00560                 int qIndex = pq%nCumulMixHar_qMax;
00561 
00562 
00563 
00564         histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] =
00565       histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq]->GetBinContent(1)/
00566       histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
00567 
00568 
00569         CpqMix[pIndex][qIndex]=mAvMult*(pow(histFull[k].histFullHar[j].mCumulIntegG0Mix[pq], 1./mAvMult) -1.);
00570 
00571 
00572           }
00573 
00574           double CpxMix[nCumulMixHar_pMax];
00575           double CpyMix[nCumulMixHar_pMax];
00576 
00577 
00578           for (int p =0; p<nCumulMixHar_pMax; p++){
00579             CpxMix[p] = (1./(4.*r0Mix))*(CpqMix[p][0] - CpqMix[p][2]);
00580             CpyMix[p] = (1./(4.*r0Mix))*(CpqMix[p][3] - CpqMix[p][1]);
00581 
00582           }
00583 
00584 
00585           cumulIntegMix[j] = (1./(4.*r0Mix*r0Mix))*(
00586                         CpxMix[0]-CpyMix[1]-CpxMix[2]+CpyMix[3]
00587                         +CpxMix[4]-CpyMix[5]-CpxMix[6]+CpyMix[7]);
00588 
00589           //use v2{4}
00590           double  tempMeanV = pow(-1.*cumulInteg2[1]*mAveMeanWgtSqr*mAveMeanWgtSqr,1./4.); // <wgt*v2>
00591           //use v2{2}
00592           //  double  tempMeanV = pow(cumulInteg1[1]*mAveMeanWgtSqr,1./2.); // <wgt*v2>
00593 
00594           cout<<"mAveMeanWgtSqr "<<mAveMeanWgtSqr<<endl;
00595           cout<<" cmumulInteg2[1] "<<cumulInteg2[1]<<endl;
00596           cout<<" tempMeanV "<<tempMeanV<<endl;
00597           double  tempMeanVMixSq = cumulIntegMix[j]/tempMeanV; //(24) 
00598  
00599 
00602             //      double chi1 = sqrt(tempMeanVMixSq)*sqrt(mAvMult); //(36)
00603             //            double chi2 = tempMeanV*sqrt(mAvMult);
00604             //            double deltaC3Sqr = (2. + 4.*chi1*chi1 + 2.*chi2*chi2 + 4.*chi1*chi1*chi2*chi2+chi1*chi1*chi1*chi1)/(2.*mAvMult*mAvMult*mAvMult*(histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1)));
00605             //
00606             //            cout<<" C3              is "<<cumulIntegMix[j]<<endl;
00607             //            cout<<" delta C3 square is "<<deltaC3Sqr<<endl;
00609  
00610 
00611         histTitle = new TString("Flow_CumulMix_v2D_Sel");
00612         histTitle->Append(*countSels);
00613         histTitle->Append("_Har");
00614         histTitle->Append(*countHars);
00615         histFull[k].histFullHar[j].mHistMix_v2D = 
00616           new TH2D(*(histFull[k].histFullHar[j].mHistCumulMix2D->
00617                      ProjectionXY(histTitle->Data(),"e")));
00618         histFull[k].histFullHar[j].mHistMix_v2D->SetTitle(histTitle->Data());
00619         histFull[k].histFullHar[j].mHistMix_v2D->SetXTitle((char*)xLabel.Data());
00620         histFull[k].histFullHar[j].mHistMix_v2D->SetYTitle("Pt (GeV)");
00621         histFull[k].histFullHar[j].mHistMix_v2D->SetZTitle("v (%)");
00622         histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./profScale);
00623 
00624         delete histTitle;
00625 
00626 
00627 
00628         histTitle = new TString("Flow_CumulMix_vEta_Sel");
00629         histTitle->Append(*countSels);
00630         histTitle->Append("_Har");
00631         histTitle->Append(*countHars);
00632         histFull[k].histFullHar[j].mHistMix_vEta = 
00633           new  TH1D(*(histFull[k].histFullHar[j].mHistCumulMixEta->
00634                       ProjectionX(histTitle->Data(),"e")));
00635         histFull[k].histFullHar[j].mHistMix_vEta->SetTitle(histTitle->Data());
00636         histFull[k].histFullHar[j].mHistMix_vEta->SetXTitle((char*)xLabel.Data());
00637         histFull[k].histFullHar[j].mHistMix_vEta->SetYTitle("v (%)");
00638         histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./profScale);
00639 
00640 
00641         delete histTitle;
00642         
00643 
00644 
00645 
00646 
00647         histTitle = new TString("Flow_CumulMix_vPt_Sel");
00648         histTitle->Append(*countSels);
00649         histTitle->Append("_Har");
00650         histTitle->Append(*countHars);
00651 
00652         //                  histTitle->Append("_TPC");
00653 
00654           histTitle2 = new TString("Flow_CumulMix2D_Sel");
00655           histTitle2->Append(*countSels);
00656           histTitle2->Append("_Har");
00657           histTitle2->Append(*countHars);
00658 
00659 
00660 
00661 
00662         histFull[k].histFullHar[j].mHistMix_vPt =
00663           //myProjectionY((TProfile2D *)f->Get(histTitle2->Data()) , 0., 1. , 1, ((j==0)?1:0) );
00664                  new  TH1D(*(histFull[k].histFullHar[j].mHistCumulMixPt->ProjectionX(histTitle->Data(),"e")));
00665 
00666         histFull[k].histFullHar[j].mHistMix_vPt->SetTitle(histTitle->Data());
00667         histFull[k].histFullHar[j].mHistMix_vPt->SetXTitle("Pt (GeV)");
00668         histFull[k].histFullHar[j].mHistMix_vPt->SetYTitle("v (%)");
00669         histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./profScale);
00670 
00671         delete histTitle;
00672         delete histTitle2;      
00673 
00674 
00675 
00676           if (tempMeanVMixSq>0.){
00677     meanIntegVMix[j] = sqrt(tempMeanVMixSq);//(24) <w1v1>
00678 
00679         histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent)); 
00680         histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
00681         histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent)); 
00682         histFull[k].mHistMix_v->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent)); 
00683        } else {
00684          histFull[k].histFullHar[j].mHistMix_v2D->Reset();
00685          histFull[k].histFullHar[j].mHistMix_vEta->Reset();
00686          histFull[k].histFullHar[j].mHistMix_vPt->Reset();
00687          histFull[k].mHistMix_v->Reset();
00688        cout<<"### <wgt*v1>**2 = "<<tempMeanVMixSq<<" < 0. 3-part mixed har ana. failed "<<endl;
00689        }
00690 
00691 
00692 
00693     } //end of Har loop for Mixed har.
00694 
00695 
00696 
00697 
00698     if (m_M==1) {
00699       
00700       TH1D* histOfMeanIntegV = new TH1D(*(histFull[k].mHist_v[0]));
00701       histOfMeanIntegV->Reset();
00702       
00703       TH1D* histOfMeanIntegV3 = new TH1D(*(histFull[k].mHist_v[1]));
00704       histOfMeanIntegV3->Reset();
00705       
00706       for (int j = 1; j < nHars+1; j++) {
00707         histOfMeanIntegV->SetBinContent(j, 1./(meanIntegV[j-1]*perCent));
00708         histOfMeanIntegV->SetBinError(j,0.);
00709         histOfMeanIntegV3->SetBinContent(j, -1./(meanIntegV3[j-1]*perCent));
00710         histOfMeanIntegV3->SetBinError(j,0.);
00711       }
00712       histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV);
00713       histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV3);
00714       
00715       //force to be zero if the value is "nan"
00716       for (int ord = 0; ord < nCumulDiffOrders; ord++){
00717         for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
00718           if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
00719                  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ) {
00720             histFull[k].mHist_v[ord]->SetBinContent(j,0.);
00721             histFull[k].mHist_v[ord]->SetBinError(j,0.);
00722           }
00723         }
00724       }
00725 
00726       for (int j = 1; j < nHars+1; j++) {
00727         if (histFull[k].mHist_v[0]->GetBinContent(j)< FLT_MAX &&
00728             histFull[k].mHist_v[0]->GetBinContent(j)> -1.*FLT_MAX) 
00729         cout << "##### 2-part v" << j << " = (" 
00730              << histFull[k].mHist_v[0]->GetBinContent(j) 
00731              <<" +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<" )"<<endl;
00732         if (histFull[k].mHist_v[1]->GetBinContent(j)< FLT_MAX &&
00733             histFull[k].mHist_v[1]->GetBinContent(j)> -1.*FLT_MAX) 
00734         cout << "##### 4-part v" << j << " = (" 
00735              << histFull[k].mHist_v[1]->GetBinContent(j) 
00736              <<" +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<" )"<<endl;
00737       }
00738       
00739       delete histOfMeanIntegV; 
00740       delete histOfMeanIntegV3;
00741       
00742     } else if (m_M==2) {
00743       
00744       TH1D* histOfMeanIntegV2 = new TH1D(*(histFull[k].mHist_v[0]));
00745       histOfMeanIntegV2->Reset();
00746       
00747       TH1D* histOfMeanIntegV4 = new TH1D(*(histFull[k].mHist_v[1]));
00748       histOfMeanIntegV4->Reset();
00749       
00750       for (int j = 1; j < nHars+1; j++) {
00751         histOfMeanIntegV2->SetBinContent(j, 1./(meanIntegV2[j-1]*perCent));
00752         histOfMeanIntegV2->SetBinError(j,0.);
00753         histOfMeanIntegV4->SetBinContent(j, -0.5/(meanIntegV4[j-1]*perCent));
00754         histOfMeanIntegV4->SetBinError(j,0.);
00755       }
00756       histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV2);
00757       histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV4);
00758 
00759       //force to be zero if the value is "nan"
00760       for (int ord = 0; ord < nCumulDiffOrders; ord++){
00761         for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
00762           if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
00763                  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ){
00764             histFull[k].mHist_v[ord]->SetBinContent(j,0.);
00765             histFull[k].mHist_v[ord]->SetBinError(j,0.);
00766           }
00767         }
00768       }
00769 
00770       for (int j = 1; j < nHars+1; j++) {
00771 
00772         if (histFull[k].mHist_v[0]->GetBinContent(j)< FLT_MAX &&
00773             histFull[k].mHist_v[0]->GetBinContent(j)> -1.*FLT_MAX) 
00774         cout << "##### 2-part v" << j << " = (" 
00775              << histFull[k].mHist_v[0]->GetBinContent(j) 
00776              <<") +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<endl;
00777         if (histFull[k].mHist_v[1]->GetBinContent(j)< FLT_MAX &&
00778             histFull[k].mHist_v[1]->GetBinContent(j)> -1.*FLT_MAX) 
00779         cout << "##### 4-part v" << j << " = (" 
00780              << histFull[k].mHist_v[1]->GetBinContent(j) 
00781              <<") +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<endl;
00782       }
00783       
00784       delete histOfMeanIntegV2; 
00785       delete histOfMeanIntegV4;
00786     }
00787 
00788   }
00789 
00790 
00791 
00792   // ============== write out ====================
00793 
00794   f->cd();
00795 
00796   for (int k = 0; k < nSels; k++) {
00797 
00798    for (int ord = 0; ord < nCumulDiffOrders; ord++) {
00799      histFull[k].mHist_v[ord]->Write(histFull[k].mHist_v[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00800    }
00801 
00802      histFull[k].mHistMix_v->Write(histFull[k].mHistMix_v->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00803 
00804 
00805     for (int j = 0; j < nHars; j++) {
00806       cout<<" writting ........................."<<endl;
00807       for (int ord = 0; ord < nCumulDiffOrders; ord++){
00808         histFull[k].histFullHar[j].mHist_v2D[ord]->Write(histFull[k].histFullHar[j].mHist_v2D[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00809         histFull[k].histFullHar[j].mHist_vEta[ord]->Write(histFull[k].histFullHar[j].mHist_vEta[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00810         histFull[k].histFullHar[j].mHist_vPt[ord]->Write(histFull[k].histFullHar[j].mHist_vPt[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00811       }
00812 
00813 
00814       if (j==0){
00815         cout<<" j "<<j<<" k "<<k<<endl;   
00816         cout<<"  histFull[k].histFullHar[j].mHistMix_v2D "<<histFull[k].histFullHar[j].mHistMix_v2D->GetTitle()<<endl;
00817 
00818       (histFull[k].histFullHar[j].mHistMix_v2D)->Write((histFull[k].histFullHar[j].mHistMix_v2D)->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00819       (histFull[k].histFullHar[j].mHistMix_vEta)->Write((histFull[k].histFullHar[j].mHistMix_vEta)->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00820       (histFull[k].histFullHar[j].mHistMix_vPt)->Write((histFull[k].histFullHar[j].mHistMix_vPt)->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
00821 
00822       }
00823 
00824    }
00825  }
00826   //=============================================
00827 
00828   f->Close();
00829 
00830 }