PeterParfenov лет назад: 4
Сommit
69012d72c9
9 измененных файлов с 4588 добавлено и 0 удалено
  1. 853 0
      FlowANA.C
  2. 1168 0
      FlowANA.h
  3. 804 0
      FlowANA_test.C
  4. 1168 0
      FlowANA_test.h
  5. 200 0
      ReadFlow.C
  6. 87 0
      ReadRes.C
  7. 28 0
      main_proc.C
  8. 28 0
      main_proc_test.C
  9. 252 0
      res2.C

+ 853 - 0
FlowANA.C

@@ -0,0 +1,853 @@
+#define FlowANA_cxx
+#include "FlowANA.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TProfile.h"
+#include "TFile.h"
+#include "TNtuple.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TRandom3.h"
+#include "TVector3.h"
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <sstream>
+#include <string>
+#include <cmath>
+
+
+
+using namespace std;
+//List of histograms and Ntuples....
+
+static const int neta = 9;
+static const int ndet = 4;
+static const int ncent = 6;
+static const int nth = 3;
+static const int npid = 4;
+
+
+
+
+static const int npt = 11; // 0.5 - 3.6 GeV/c - number of pT bins
+static const double bin_w[11]={0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.8,2.3,2.8,4.0};
+
+static const float maxpt = 4.0; // max pt
+static const float minpt = 0.2; // min pt
+
+
+TH1F *hpt[ncent][npt][npid];
+TH1F *hv2[ndet][ncent][npt][npid];
+TH1F *hv22[ndet][npt][npid];
+
+
+
+TH1F *H_Qw[neta];
+TH1F *H_EP[nth][neta];
+TH1F *H_Qv[nth][neta];
+
+
+TH1F *HRes[nth][ndet][ncent];
+
+// for centrality determination
+TH1F *href1; // PHENIX BBC
+TH1F *href2; // PHENIX RXN
+TH1F *href3; // STAR TPC
+TH1F *href4; // STAR ntracks>2
+TH1F *href5; // STAR ntracks>4
+
+
+TH1F *hbimp; // impact parameter
+TH1F *hbimp2; // impact parameter
+TH1F *hbimp3; // impact parameter
+TH1F *hbimp4; // impact parameter
+TH1F *hbimp5; // impact parameter
+
+TH1F *h2pt;
+TH1F *h2eta;
+TH1F *h2phi;
+TH1F *h2phis;
+TH1F *h2phirp;
+
+
+
+
+static const double MYPI=3.141592654;
+
+TFile *d_outfile;
+
+
+void FlowANA::Loop()
+{
+
+   if (fChain == 0) return;
+
+   Int_t nentries = Int_t(fChain->GetEntries());
+
+   Int_t nbytes = 0, nb = 0;
+ 
+
+   for (Int_t jentry=0; jentry<nentries;jentry++) {
+
+    
+      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
+      nb = fChain->GetEntry(jentry);   nbytes += nb;
+      
+      if(jentry%1000==0) cout << "Event [" << jentry << "/" << nentries << "]" <<endl;// event counter
+     
+      ana_event(jentry,ientry);
+   }
+}
+
+void FlowANA::ana_init(const char *outfile) {
+
+book_hist(outfile);
+ gRandom->SetSeed( (unsigned int)time(NULL) ); 
+}
+
+
+
+
+void FlowANA::loop_a_file(const char *file) {
+  TFile *treefile = TFile::Open(file);
+  TTree *tree = (TTree*)treefile->Get("mctree"); // put the name of the TTree
+  if(tree == 0) {
+    cout <<"htree is not found in "<<file<<endl;
+    treefile->Close();
+    return;
+  }
+  cout << file <<" is opened"<<endl;
+  Init(tree);
+
+  
+
+  Loop();
+  treefile->Close();
+  cout <<"one file processed"<<endl;
+}
+
+void FlowANA::ana_end() {
+d_outfile->cd();
+
+for( int ieta=0; ieta<neta; ieta++ ){
+		H_Qw[ieta]->Write();
+	}
+
+
+ 	for( int ith=0; ith<3; ith++ ){
+		for( int ieta=0; ieta<neta; ieta++ ){
+			H_EP[ith][ieta]->Write();
+			H_Qv[ith][ieta]->Write();
+		}
+	}
+
+
+	for( int ith=0; ith<3; ith++ ){
+		for( int idet=0; idet<ndet; idet++ ){
+                         for( int icent=0; icent<ncent; icent++ ){
+		  
+			HRes[ith][idet][icent]->Write();
+		    
+		}
+	}
+
+}
+	
+
+for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+			            hpt[icent][ipt][id]->Write();
+		    
+		         }
+	            }
+
+             }
+
+
+for( int idet=0; idet<ndet; idet++ ){
+               for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+			            hv2[idet][icent][ipt][id]->Write();
+		    
+		         }
+	            }
+
+             }
+
+}
+
+
+for( int idet=0; idet<ndet; idet++ ){
+               
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+			            hv22[idet][ipt][id]->Write();
+		    
+		         
+	            }
+
+             }
+
+}
+
+
+	
+h2pt->Write();
+h2eta->Write();
+h2phis->Write();
+h2phi->Write();
+h2phirp->Write();
+
+ 
+href1->Write();
+href2->Write();
+href3->Write();
+href4->Write();
+href5->Write();
+
+
+ 
+hbimp->Write();
+hbimp2->Write();
+hbimp3->Write(); 
+hbimp4->Write();
+hbimp5->Write(); 
+
+ d_outfile->Close();
+
+
+
+} // end of ana_end
+
+void FlowANA::book_hist(const char *outfile) {
+  d_outfile = new TFile(outfile,"recreate");
+
+
+
+h2eta  = new TH1F("h2eta","eta distribution",200,-7,7);
+h2pt  = new TH1F("h2pt","pt distribution",200,0,10.0);
+
+
+ 
+href1  = new TH1F("href1","Centrality by BBC (PHENIX)",1000,0,1000);
+href2  = new TH1F("href2","Centrality by BBC (PHENIX)",1000,0,1000);
+href3  = new TH1F("href3","Centrality by TPC (STAR)",1000,0,1000);
+href4  = new TH1F("href4","Centrality by TPC (STAR), ntracks >1",1000,0,1000);
+href5  = new TH1F("href5","Centrality by TPC (STAR), ntracks >2",1000,0,1000);
+
+hbimp2  = new TH1F("hbimp2","Imact Parameter",180,0,18);
+hbimp3  = new TH1F("hbimp3","Imact Parameter",180,0,18);
+hbimp4  = new TH1F("hbimp4","Imact Parameter",180,0,18);
+hbimp5  = new TH1F("hbimp5","Imact Parameter",180,0,18);
+
+ 
+hbimp  = new TH1F("hbimp","Imact Parameter",180,0,18);
+
+h2phi = new TH1F("h2phi","azim. angle", 100, -TMath::Pi(), TMath::Pi()+2.*TMath::Pi()/100 );
+h2phis = new TH1F("h2phis","azim. angle (|eta| <0.5)", 100, -TMath::Pi(), TMath::Pi()+2.*TMath::Pi()/100 );
+h2phirp = new TH1F("h2phirp","azim. angle", 100, -TMath::Pi(), TMath::Pi()+2.*TMath::Pi()/100 );
+
+
+for( int ieta=0; ieta<neta; ieta++ ) H_Qw[ieta] = new TH1F( Form("H_Qw_%d",ieta), Form("H_Qw_%d",ieta), 100, 0, 1000 );
+  
+
+for( int ith=0; ith<3; ith++ ){
+     for( int ieta=0; ieta<neta; ieta++ ){
+  H_EP[ith][ieta] = new TH1F( Form("H_EP_%d_%d",ith,ieta), Form("H_EP_%d_%d",ith,ieta), 100, -TMath::Pi()/(ith+2.)-0.1, TMath::Pi()/(ith+2.)+0.10 );
+ H_Qv[ith][ieta] = new TH1F( Form("H_Qv_%d_%d",ith,ieta), Form("H_Qv_%d_%d",ith,ieta), 100, 0, 10 );
+		}
+	}
+
+
+
+	for( int ith=0; ith<3; ith++ ){
+		for( int idet=0; idet<ndet; idet++ ){
+                         for( int icent=0; icent<ncent; icent++ ){
+   HRes[ith][idet][icent] = new TH1F( Form("HRes_%d_%d_%d",ith,idet,icent), Form("HRes_%d_%d_%d",ith,idet,icent), 100, -10, 10 );	       
+		    
+		}
+	}
+
+}
+
+
+for( int idet=0; idet<ndet; idet++ ){
+               for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+hv2[idet][icent][ipt][id] = new TH1F( Form("hv2_%d_%d_%d_%d",idet,icent,ipt,id), Form("hv2_%d_%d_%d_%d",idet,icent,ipt,id), 400, -10, 10 );			    
+		    
+		         }
+	            }
+
+             }
+
+}
+
+
+
+for( int idet=0; idet<ndet; idet++ ){
+             
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+hv22[idet][ipt][id] = new TH1F( Form("hv22_%d_%d_%d",idet,ipt,id), Form("hv22_%d_%d_%d",idet,ipt,id), 400, -10, 10 );			    
+		    
+		       
+	            }
+
+             }
+
+}
+
+
+
+               for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+hpt[icent][ipt][id] = new TH1F( Form("hpt_%d_%d_%d",icent,ipt,id), Form("hpt_%d_%d_%d",icent,ipt,id), 400, 0, 10 );			    
+		    
+		         }
+	            }
+
+             }
+
+
+
+
+	
+
+}
+
+
+
+
+
+
+
+// Analysis for each event !!!!  
+
+void FlowANA::ana_event(int jentry, int ientry) { 
+
+
+  float phiRP = gRandom->Uniform(0, 2.*TMath::Pi());
+
+  // float phiRP = gRandom->Uniform(-1.0*TMath::Pi(),TMath::Pi());
+ h2phirp->Fill(phiRP);
+
+  // centrality cut and vertex +/- 30 cm cut
+  /*
+
+ if(cent>0&&cent<=80){
+
+
+       if(centrality<=5)        mycent=0;
+       else if(centrality<=10)  mycent=1;
+       else if(centrality<=15)  mycent=2;
+       else if(centrality<=20)  mycent=3;
+       else if(centrality<=25)  mycent=4;
+       else if(centrality<=30)  mycent=5;
+       else if(centrality<=35)  mycent=6;
+       else if(centrality<=40)  mycent=7;
+       else if(centrality<=45)  mycent=8;
+       else if(centrality<=50)  mycent=9;
+       else if(centrality<=55)  mycent=10;
+
+      else   mycent=11;
+
+  */
+
+ //if(ientry%100000==0) cout << ientry <<endl;// event counter
+
+
+
+	float sumQxy[3][9][2] = {{{0}}}; //[ith][eta][x,y]
+	//float sumQxyFull[3][2] = {{{0}}};
+	float multQv[9] = {0}; //[eta]
+
+
+ 
+hbimp->Fill(bimp);
+
+double refMult1 = 0;
+double refMult2 = 0;
+
+ 
+
+double Nch_L = 0;
+double Nch_R = 0;
+
+double Nch_L2 = 0;
+double Nch_R2 = 0;
+ 
+  for(int itrk=0;itrk<nh;itrk++) {  //track loop
+
+        float pt  = sqrt( TMath::Power(momx[itrk], 2.0 ) + TMath::Power(momy[itrk], 2.0 ) );
+	
+	float oldphi = atan2( momx[itrk], momy[itrk] );
+	float phi=oldphi;
+	
+	float the = atan2( pt, momz[itrk] );//atan2(pt/pz)
+	float eta = -log( tan( 0.5 * the ) );
+
+	//	phi += phiRP;
+	// float px = pt * cos(phi);
+	//	float py = pt * sin(phi);
+	//	pt  = sqrt( TMath::Power(px, 2.0 ) + TMath::Power(py, 2.0 ) );
+
+            h2pt->Fill(pt);
+            h2eta->Fill(eta);
+            h2phi->Fill(oldphi);
+            h2phis->Fill(phi);
+
+	    
+        if( pt>0.1 && fabs(eta)<0.5 ) refMult1++;
+	if( pt>0.0 && fabs(eta)<0.5 ) refMult2++;
+        if(eta >=  3.1 && eta <=  4.0) Nch_R++;
+        if(eta >= -4.0 && eta <= -3.1) Nch_L++;
+
+       
+        if( pt<0.15  || pt>2. ) continue;
+
+        int fEta = -1;
+
+	// TPC plane
+	if( eta>-1 && eta<-0.1 ) fEta = 0; // TPC East
+	if( eta>0.1 && eta<1 )   fEta = 1; // TPC West
+	if( fabs(eta)<0.1 )      fEta = 2; // TPC Mid
+
+      	// RXN plane
+       	if( eta>-3.0 && eta<-1.0 )   fEta = 3; // RXN East
+       	if( eta>1.0 && eta<3.0 )     fEta = 4; // RXN West
+
+	// BBC plane 
+	if( eta>-5 && eta<-3 )      fEta = 5; //East
+	if( eta>3.0 && eta<5  )     fEta = 6; //West
+
+	// FHCal plane
+	if( eta>-5 && eta<-2 )      fEta = 7; //East
+	if( eta>2.0 && eta<5  )     fEta = 8; //West
+
+	// if( fabs(eta)>1.1 && fabs(eta)<2.9 )     fEta = 7; // RXN combined
+        //if( fabs(eta)>3.0 && fabs(eta)<5.0 )     fEta = 8; // BBC combined
+	
+
+
+           if( fEta>-1 ){
+			for( int ith=0; ith<3; ith++ ){
+				if (fEta < 7){
+					sumQxy[ith][fEta][0] += pt * cos( (ith+2.0) * phi );
+					sumQxy[ith][fEta][1] += pt * sin( (ith+2.0) * phi );
+				}else{
+					sumQxy[ith][fEta][0] += pt * cos( (ith+1.0) * phi );
+					sumQxy[ith][fEta][1] += pt * sin( (ith+1.0) * phi );
+				}
+				//if(fEta == 7 || fEta == 8){
+				//	sumQxyFull[ith][0] += pt * cos( (ith+1.0) * phi );
+				//	sumQxyFull[ith][1] += pt * sin( (ith+1.0) * phi );
+				//}
+			}
+			multQv[fEta]++;
+	   }// end of eta selection
+      
+
+   }// end of track loop
+
+
+  
+      
+      if(Nch_L >= 1 && Nch_R >= 1){
+          href1 -> Fill(Nch_L + Nch_R);
+      }
+
+      if(Nch_L >= 2 && Nch_R >= 2){
+          href2 -> Fill(Nch_L + Nch_R);
+	  hbimp2->Fill(bimp);
+      }
+      if(refMult1>0){
+          href3 -> Fill(refMult1);
+	  hbimp3->Fill(bimp);
+      }
+      
+
+      if(refMult2>0){
+          href4 -> Fill(refMult1);
+	  hbimp4->Fill(bimp);
+      }
+
+      if(refMult1>2){
+          href5 -> Fill(refMult1);
+	  hbimp5->Fill(bimp);
+      }
+
+      float sumLR=Nch_L + Nch_R;
+  
+      // int fCent   = GetCentrality10_RefMult( refMult1 );// STAR def
+      //if( fCent<0 ) cout << fCent << endl;
+
+      // int fCent   =  GetCentrality10_RefMultPHENIX(sumLR);
+
+      int fCent   = GetCentrality10_Bimp(bimp);
+      //int fCent   =  GetCentrality10_BimpExp(bimp);
+      
+  	float fEP[3][9]; //[ith][eta]
+	float fQv[3][9]; //[ith][eta]
+	for( int ith=0; ith<3; ith++ ){ // flow harmonic loop
+	  for( int ieta=0; ieta<9; ieta++ ){ // ep detector gap 
+	    if( multQv[ieta]>5 ){ // multiplicity > 5
+			if( ieta<7 )
+			{
+				fEP[ith][ieta] = atan2( sumQxy[ith][ieta][1], sumQxy[ith][ieta][0] ) / ( ith + 2.0 );
+				fEP[ith][ieta] = atan2( sin( (ith+2.0)*fEP[ith][ieta] ), cos( (ith+2.0)*fEP[ith][ieta] ) );
+				fEP[ith][ieta] /= ( ith + 2.0 );
+			}else{
+				fEP[ith][ieta] = atan2( sumQxy[ith][ieta][1], sumQxy[ith][ieta][0] ) / ( ith + 1.0 );
+				fEP[ith][ieta] = atan2( sin( (ith+1.0)*fEP[ith][ieta] ), cos( (ith+1.0)*fEP[ith][ieta] ) );
+				fEP[ith][ieta] /= ( ith + 1.0 );
+			}
+ fQv[ith][ieta] = sqrt( TMath::Power( sumQxy[ith][ieta][0],2.0)+TMath::Power( sumQxy[ith][ieta][1], 2.0))/sqrt( multQv[ieta]);
+ H_Qw[ieta]->Fill( multQv[ieta] );
+		}else{
+			fEP[ith][ieta] = -9999;
+			fQv[ith][ieta] = -9999;
+	    }// end of mult cut selection
+	  } // end of loop on EP detectors
+	} // end of flow harmonic loop
+
+
+
+	for( int ith=0; ith<3; ith++ ){ // harmonic loop
+	  for( int ieta=0; ieta<9; ieta++ ){// eta EP detector loop
+	    if( fEP[ith][ieta]>-9000 ){ // EP reconstructed 
+				H_EP[ith][ieta]->Fill( fEP[ith][ieta] );
+				H_Qv[ith][ieta]->Fill( fQv[ith][ieta] );
+	    }// end of EP reconstructed
+	  }// end of eta loop
+	}// end of harm loop
+
+	
+		//Resolution
+	for( int ith=0; ith<3; ith++ ){
+		for( int icb=0; icb<4; icb++ ){
+			double psi1, psi2, fq1, fq2;
+       
+	       if ( icb==0 ){ psi1 = fEP[ith][0]; psi2 = fEP[ith][1]; fq1 = fQv[ith][0]; fq2 = fQv[ith][1]; } // TPC.E-TPC.W
+	       else if( icb==1 ){ psi1 = fEP[ith][3]; psi2 = fEP[ith][4]; fq1 = fQv[ith][3];  fq2 = fQv[ith][4]; } // RXN.E-RXN.W
+	       else if( icb==2 ){ psi1 = fEP[ith][5]; psi2 = fEP[ith][6]; fq1 = fQv[ith][5]; fq2 = fQv[ith][6]; } // BBC.E-BBC.W
+	       else { psi1 = fEP[ith][7]; psi2 = fEP[ith][8]; fq1 = fQv[ith][7]; fq2 = fQv[ith][8]; } // FHCal.E-FHCal.W
+		       
+
+		        
+		
+		
+
+			if( psi1<-9000 || psi2<-9000 ) continue;
+			if( fq1<0 || fq2<0 ) continue;
+
+
+
+                        double dPsi = ( ith + 2. ) * ( psi1 - psi2 );
+		   	dPsi = atan2( sin(dPsi), cos(dPsi) );
+			if(fCent>-1&&fCent<6){
+			
+			  HRes[ith][icb][fCent]->Fill(cos(dPsi) );
+
+			}
+
+		}
+	}
+
+
+
+// refmult star
+//float res2tpc[6]={0.503463,0.71591,0.749962,0.708934,0.61689,0.475386};
+//float res2rxn[6]={0.547883,0.791309,0.824213,0.793991,0.709638,0.561604};
+//float res2bbc[6]={0.252904,0.377894,0.401809,0.366931,0.300508,0.223784};
+
+// phenix ala cent
+
+// float res2tpc[6]={0.482228,0.704061,0.749667,0.7259,0.655706,0.541163};
+//float res2rxn[6]={ 0.526791,0.780304,0.824815,0.806745,0.746894,0.634004};
+//float res2bbc[6]={0.241286,0.375359,0.40532,0.384438,0.326188,0.250811};
+
+//float res2tpc[6]={0.493507,0.716751,0.748541,0.709883,0.623788,0.491104};
+
+//float res2rxn[6]={0.540663,0.793812,0.82354,0.794412,0.714603,0.579696};
+//float res2bbc[6]={0.248559,0.383837,0.40206,0.373202,0.299496,0.2273};
+//float res2fhcal[6]={0.38477,0.577671,0.601972,0.555919,0.46528,0.339993};
+
+// 11.5 gev
+//
+float res2tpc[6] = {0.323761,0.434068,0.433521,0.383028,0.307532,0.250556};
+float res2rxn[6] = {0.212259,0.333732,0.340171,0.303697,0.247778,0.206723};
+float res2bbc[6] = {1,1,1,1,1,1};
+float res2fhcal[6] = {0.219457,0.37629,0.417876,0.39807,0.340169,0.271286};
+
+// 7.7 gev
+//
+//float res2tpc[6] = {0.279143,0.404904,0.410895,0.364538,0.298188,0.249491};
+//float res2rxn[6] = {0.171143,0.278782,0.277882,0.234889,0.180159,0.142501};
+//float res2bbc[6] = {0,0,0,0,0,0};
+//float res2fhcal[6] = {0.261058,0.463961,0.520235,0.504925,0.444674,0.355923};
+
+if(fCent>=0&&fCent<6){
+
+ 
+
+  
+  for(int itrk=0;itrk<nh;itrk++) {  //track loop
+
+        float pt  = sqrt( TMath::Power(momx[itrk], 2.0 ) + TMath::Power(momy[itrk], 2.0 ) );
+	
+	float oldphi = atan2( momx[itrk], momy[itrk] );
+	float phi=oldphi;
+	
+	float the = atan2( pt, momz[itrk] );//atan2(pt/pz)
+	float eta = -log( tan( 0.5 * the ) );
+
+        // phi += phiRP;
+	//  float px = pt * cos(phi);
+	//	float py = pt * sin(phi);
+	//	pt  = sqrt( TMath::Power(px, 2.0 ) + TMath::Power(py, 2.0 ) );
+
+
+	
+
+        if( pt<0.15  || pt>4.0 ) continue;
+
+
+	
+       int ipt = 0;
+ 	for(int i=0; i<npt-1;i++){
+	if(pt>=bin_w[i] && pt<bin_w[i+1]) ipt = i;
+	}// end of ipt loop
+
+        float v2tpc=-999.0;
+	float v2rxn=-999.0;
+	float v2bbc=-999.0;
+	float v2fhcal=-999.0;
+
+	if(eta>0&&eta<1.0){
+         v2tpc = cos(2.0 * (phi-fEP[0][0]) )/res2tpc[fCent];
+         v2rxn = cos(2.0 * (phi-fEP[0][3]) )/res2rxn[fCent];
+         v2bbc = cos(2.0 * (phi-fEP[0][5]) )/res2bbc[fCent];
+         v2fhcal = cos(2.0 * (phi-fEP[0][7]) )/res2fhcal[fCent];
+ 
+	}
+
+        if(eta<0&&eta>-1.0){
+         v2tpc = cos(2.0 * (phi-fEP[0][1]) )/res2tpc[fCent];
+         v2rxn = cos(2.0 * (phi-fEP[0][4]) )/res2rxn[fCent];
+         v2bbc = cos(2.0 * (phi-fEP[0][6]) )/res2bbc[fCent];
+         v2fhcal = cos(2.0 * (phi-fEP[0][8]) )/res2fhcal[fCent];
+ 
+	}
+
+        
+	
+
+        if(fabs(eta)<1.0){
+	  hv2[0][fCent][ipt][0]->Fill(v2tpc);
+          hv2[1][fCent][ipt][0]->Fill(v2rxn);
+          hv2[2][fCent][ipt][0]->Fill(v2bbc);
+          hv2[3][fCent][ipt][0]->Fill(v2fhcal);
+	  hpt[fCent][ipt][0]->Fill(pt);
+
+          if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][0]->Fill(v2tpc);
+          hv22[1][ipt][0]->Fill(v2rxn);
+          hv22[2][ipt][0]->Fill(v2bbc);
+          hv22[3][ipt][0]->Fill(v2fhcal);
+
+	  }
+	  
+
+	  if ( pdg[itrk]==211){
+              hv2[0][fCent][ipt][1]->Fill(v2tpc);
+              hv2[1][fCent][ipt][1]->Fill(v2rxn);
+              hv2[2][fCent][ipt][1]->Fill(v2bbc);
+			  hv2[3][fCent][ipt][1]->Fill(v2fhcal);
+	      hpt[fCent][ipt][1]->Fill(pt);
+
+
+          if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][1]->Fill(v2tpc);
+          hv22[1][ipt][1]->Fill(v2rxn);
+          hv22[2][ipt][1]->Fill(v2bbc);
+          hv22[3][ipt][1]->Fill(v2fhcal);
+
+	  }
+
+
+	      
+
+          }// end of pion selection
+
+
+	  if ( pdg[itrk]==321){
+              hv2[0][fCent][ipt][2]->Fill(v2tpc);
+              hv2[1][fCent][ipt][2]->Fill(v2rxn);
+              hv2[2][fCent][ipt][2]->Fill(v2bbc);
+			  hv2[3][fCent][ipt][2]->Fill(v2fhcal);
+	      hpt[fCent][ipt][2]->Fill(pt);
+
+          if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][2]->Fill(v2tpc);
+          hv22[1][ipt][2]->Fill(v2rxn);
+          hv22[2][ipt][2]->Fill(v2bbc);
+          hv22[3][ipt][2]->Fill(v2fhcal);
+
+	  }
+	      
+
+          }// end of kaon selection
+
+
+	   if ( pdg[itrk]==2212){
+              hv2[0][fCent][ipt][3]->Fill(v2tpc);
+              hv2[1][fCent][ipt][3]->Fill(v2rxn);
+              hv2[2][fCent][ipt][3]->Fill(v2bbc);
+			  hv2[3][fCent][ipt][3]->Fill(v2fhcal);
+	      hpt[fCent][ipt][3]->Fill(pt);
+
+           if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][3]->Fill(v2tpc);
+          hv22[1][ipt][3]->Fill(v2rxn);
+          hv22[2][ipt][3]->Fill(v2bbc);
+          hv22[3][ipt][3]->Fill(v2fhcal);
+
+	  }
+	      
+
+          }// end of proton selection
+
+	   
+	}
+
+        
+
+
+
+
+
+	
+
+
+  }// end of the track loop
+
+ }// end of centrality selection 
+
+
+
+
+
+
+
+
+
+	
+     
+ }
+
+
+
+// Urqmd 7.7 GeV
+
+int FlowANA::GetCentrality10_RefMult( double refMult ){
+
+	int fcent;
+	if     ( refMult>=268 ) fcent = 0; // 0-10% 
+	else if( refMult>=185 ) fcent = 1; //10-20% 
+	else if( refMult>=126 ) fcent = 2; //20-30% 
+	else if( refMult>=83 ) fcent = 3; //30-40% 
+	else if( refMult>=52 ) fcent = 4; //40-50% 
+	else if( refMult>= 21 ) fcent = 5; //50-60% 
+	else if( refMult>= 17 ) fcent = 6; //60-70% 
+	else if( refMult>= 9 ) fcent = 7; //70-80% 
+	else                    fcent =-1;
+
+	return fcent;
+}
+
+
+int FlowANA::GetCentrality10_RefMultPHENIX( double refMult ){
+
+	int fcent;
+	if     ( refMult>=516 ) fcent = 0; // 0-10% 
+	else if( refMult>=382 ) fcent = 1; //10-20% 
+	else if( refMult>=280 ) fcent = 2; //20-30% 
+	else if( refMult>=199 ) fcent = 3; //30-40% 
+	else if( refMult>=137 ) fcent = 4; //40-50% 
+	else if( refMult>= 90 ) fcent = 5; //50-60% 
+	else if( refMult>= 55 ) fcent = 6; //60-70% 
+	else if( refMult>= 31 ) fcent = 7; //70-80% 
+	else                    fcent =-1;
+
+	return fcent;
+}
+
+
+
+
+
+
+
+int FlowANA::GetCentrality10_Bimp( float bimp ){
+
+	int fcent;
+	if     ( bimp<4.5 ) fcent = 0; // 0-10%
+	else if( bimp<6.3 )  fcent = 1; //10-20%
+	else if( bimp<7.73 ) fcent = 2; //20-30%
+	else if( bimp<8.92 ) fcent = 3; //30-40%
+	else if( bimp<9.99)  fcent = 4; //40-50%
+	else if( bimp<10.83) fcent = 5; //50-60%
+	else if( bimp<11.78) fcent = 6; //60-70%
+	else if( bimp<12.61) fcent = 7; //70-80%
+	else                fcent =-1;
+
+	return fcent;
+}
+
+
+int FlowANA::GetCentrality10_BimpExp( float bimp ){
+
+	int fcent;
+	if     ( bimp<4.68 ) fcent = 0; // 0-10%
+	else if( bimp<6.47 ) fcent = 1; //10-20%
+	else if( bimp<7.99 ) fcent = 2; //20-30%
+	else if( bimp<9.31 ) fcent = 3; //30-40%
+	else if( bimp<10.48 ) fcent = 4; //40-50%
+	else if( bimp<11.49 ) fcent = 5; //50-60%
+	else if( bimp<12.35) fcent = 6; //60-70%
+	else if( bimp<13.0) fcent = 7; //70-80%
+	else                fcent =-1;
+
+	return fcent;
+}
+

Разница между файлами не показана из-за своего большого размера
+ 1168 - 0
FlowANA.h


+ 804 - 0
FlowANA_test.C

@@ -0,0 +1,804 @@
+#define FlowANA_cxx
+#include "FlowANA_test.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TProfile.h"
+#include "TFile.h"
+#include "TNtuple.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TRandom3.h"
+#include "TVector3.h"
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <sstream>
+#include <string>
+#include <cmath>
+
+
+
+using namespace std;
+//List of histograms and Ntuples....
+
+static const int neta = 7;
+static const int ndet = 3;
+static const int ncent = 6;
+static const int nth = 3;
+static const int npid = 4;
+
+
+
+
+static const int npt = 11; // 0.5 - 3.6 GeV/c - number of pT bins
+static const double bin_w[11]={0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.8,2.3,2.8,4.0};
+
+static const float maxpt = 4.0; // max pt
+static const float minpt = 0.2; // min pt
+
+
+TH1F *hpt[ncent][npt][npid];
+TH1F *hv2[ndet][ncent][npt][npid];
+TH1F *hv22[ndet][npt][npid];
+
+
+
+TH1F *H_Qw[neta];
+TH1F *H_EP[nth][neta];
+TH1F *H_Qv[nth][neta];
+
+
+TH1F *HRes[nth][ndet][ncent];
+
+// for centrality determination
+TH1F *href1; // PHENIX BBC
+TH1F *href2; // PHENIX RXN
+TH1F *href3; // STAR TPC
+TH1F *href4; // STAR ntracks>2
+TH1F *href5; // STAR ntracks>4
+
+
+TH1F *hbimp; // impact parameter
+TH1F *hbimp2; // impact parameter
+TH1F *hbimp3; // impact parameter
+TH1F *hbimp4; // impact parameter
+TH1F *hbimp5; // impact parameter
+
+TH1F *h2pt;
+TH1F *h2eta;
+TH1F *h2phi;
+TH1F *h2phis;
+TH1F *h2phirp;
+
+
+
+
+static const double MYPI=3.141592654;
+
+TFile *d_outfile;
+
+
+void FlowANA::Loop()
+{
+
+   if (fChain == 0) return;
+
+   Int_t nentries = Int_t(fChain->GetEntries());
+
+   Int_t nbytes = 0, nb = 0;
+ 
+
+   for (Int_t jentry=0; jentry<nentries;jentry++) {
+
+    
+      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
+      nb = fChain->GetEntry(jentry);   nbytes += nb;
+     
+      ana_event(jentry,ientry);
+   }
+}
+
+void FlowANA::ana_init(char *outfile) {
+
+book_hist(outfile);
+ gRandom->SetSeed( (unsigned int)time(NULL) ); 
+}
+
+
+
+
+void FlowANA::loop_a_file(char *file) {
+  TFile *treefile = TFile::Open(file);
+  TTree *tree = (TTree*)treefile->Get("mctree"); // put the name of the TTree
+  if(tree == 0) {
+    cout <<"htree is not found in "<<file<<endl;
+    treefile->Close();
+    return;
+  }
+  cout << file <<" is opened"<<endl;
+  Init(tree);
+
+  
+
+  Loop();
+  treefile->Close();
+  cout <<"one file processed"<<endl;
+}
+
+void FlowANA::ana_end() {
+d_outfile->cd();
+
+for( int ieta=0; ieta<neta; ieta++ ){
+		H_Qw[ieta]->Write();
+	}
+
+
+ 	for( int ith=0; ith<3; ith++ ){
+		for( int ieta=0; ieta<neta; ieta++ ){
+			H_EP[ith][ieta]->Write();
+			H_Qv[ith][ieta]->Write();
+		}
+	}
+
+
+	for( int ith=0; ith<3; ith++ ){
+		for( int idet=0; idet<ndet; idet++ ){
+                         for( int icent=0; icent<ncent; icent++ ){
+		  
+			HRes[ith][idet][icent]->Write();
+		    
+		}
+	}
+
+}
+	
+
+for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+			            hpt[icent][ipt][id]->Write();
+		    
+		         }
+	            }
+
+             }
+
+
+for( int idet=0; idet<ndet; idet++ ){
+               for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+			            hv2[idet][icent][ipt][id]->Write();
+		    
+		         }
+	            }
+
+             }
+
+}
+
+
+for( int idet=0; idet<ndet; idet++ ){
+               
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+			            hv22[idet][ipt][id]->Write();
+		    
+		         
+	            }
+
+             }
+
+}
+
+
+	
+h2pt->Write();
+h2eta->Write();
+h2phis->Write();
+h2phi->Write();
+h2phirp->Write();
+
+ 
+href1->Write();
+href2->Write();
+href3->Write();
+href4->Write();
+href5->Write();
+
+
+ 
+hbimp->Write();
+hbimp2->Write();
+hbimp3->Write(); 
+hbimp4->Write();
+hbimp5->Write(); 
+
+ d_outfile->Close();
+
+
+
+} // end of ana_end
+
+void FlowANA::book_hist(char *outfile) {
+  d_outfile = new TFile(outfile,"recreate");
+
+
+
+h2eta  = new TH1F("h2eta","eta distribution",200,-7,7);
+h2pt  = new TH1F("h2pt","pt distribution",200,0,10.0);
+
+
+ 
+href1  = new TH1F("href1","Centrality by BBC (PHENIX)",1000,0,1000);
+href2  = new TH1F("href2","Centrality by BBC (PHENIX)",1000,0,1000);
+href3  = new TH1F("href3","Centrality by TPC (STAR)",1000,0,1000);
+href4  = new TH1F("href4","Centrality by TPC (STAR), ntracks >1",1000,0,1000);
+href5  = new TH1F("href5","Centrality by TPC (STAR), ntracks >2",1000,0,1000);
+
+hbimp2  = new TH1F("hbimp2","Imact Parameter",180,0,18);
+hbimp3  = new TH1F("hbimp3","Imact Parameter",180,0,18);
+hbimp4  = new TH1F("hbimp4","Imact Parameter",180,0,18);
+hbimp5  = new TH1F("hbimp5","Imact Parameter",180,0,18);
+
+ 
+hbimp  = new TH1F("hbimp","Imact Parameter",180,0,18);
+
+h2phi = new TH1F("h2phi","azim. angle", 100, -TMath::Pi(), TMath::Pi()+2.*TMath::Pi()/100 );
+h2phis = new TH1F("h2phis","azim. angle (|eta| <0.5)", 100, -TMath::Pi(), TMath::Pi()+2.*TMath::Pi()/100 );
+h2phirp = new TH1F("h2phirp","azim. angle", 100, -TMath::Pi(), TMath::Pi()+2.*TMath::Pi()/100 );
+
+
+for( int ieta=0; ieta<neta; ieta++ ) H_Qw[ieta] = new TH1F( Form("H_Qw_%d",ieta), Form("H_Qw_%d",ieta), 100, 0, 1000 );
+  
+
+for( int ith=0; ith<3; ith++ ){
+     for( int ieta=0; ieta<neta; ieta++ ){
+  H_EP[ith][ieta] = new TH1F( Form("H_EP_%d_%d",ith,ieta), Form("H_EP_%d_%d",ith,ieta), 100, -TMath::Pi()/(ith+2.)-0.1, TMath::Pi()/(ith+2.)+0.10 );
+ H_Qv[ith][ieta] = new TH1F( Form("H_Qv_%d_%d",ith,ieta), Form("H_Qv_%d_%d",ith,ieta), 100, 0, 10 );
+		}
+	}
+
+
+
+	for( int ith=0; ith<3; ith++ ){
+		for( int idet=0; idet<ndet; idet++ ){
+                         for( int icent=0; icent<ncent; icent++ ){
+   HRes[ith][idet][icent] = new TH1F( Form("HRes_%d_%d_%d",ith,idet,icent), Form("HRes_%d_%d_%d",ith,idet,icent), 100, -10, 10 );	       
+		    
+		}
+	}
+
+}
+
+
+for( int idet=0; idet<ndet; idet++ ){
+               for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+hv2[idet][icent][ipt][id] = new TH1F( Form("hv2_%d_%d_%d_%d",idet,icent,ipt,id), Form("hv2_%d_%d_%d_%d",idet,icent,ipt,id), 400, -10, 10 );			    
+		    
+		         }
+	            }
+
+             }
+
+}
+
+
+
+for( int idet=0; idet<ndet; idet++ ){
+             
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+hv22[idet][ipt][id] = new TH1F( Form("hv22_%d_%d_%d",idet,ipt,id), Form("hv22_%d_%d_%d",idet,ipt,id), 400, -10, 10 );			    
+		    
+		       
+	            }
+
+             }
+
+}
+
+
+
+               for( int icent=0; icent<ncent; icent++ ){
+
+			   for( int ipt=0; ipt<npt-1; ipt++ ){
+
+			     for( int id=0; id<npid; id++ ){
+		  
+hpt[icent][ipt][id] = new TH1F( Form("hpt_%d_%d_%d",icent,ipt,id), Form("hpt_%d_%d_%d",icent,ipt,id), 400, 0, 10 );			    
+		    
+		         }
+	            }
+
+             }
+
+
+
+
+	
+
+}
+
+
+
+
+
+
+
+// Analysis for each event !!!!  
+
+void FlowANA::ana_event(int jentry, int ientry) { 
+
+
+  float phiRP = gRandom->Uniform(0, 2.*TMath::Pi());
+
+  // float phiRP = gRandom->Uniform(-1.0*TMath::Pi(),TMath::Pi());
+ h2phirp->Fill(phiRP);
+
+  // centrality cut and vertex +/- 30 cm cut
+  /*
+
+ if(cent>0&&cent<=80){
+
+
+       if(centrality<=5)        mycent=0;
+       else if(centrality<=10)  mycent=1;
+       else if(centrality<=15)  mycent=2;
+       else if(centrality<=20)  mycent=3;
+       else if(centrality<=25)  mycent=4;
+       else if(centrality<=30)  mycent=5;
+       else if(centrality<=35)  mycent=6;
+       else if(centrality<=40)  mycent=7;
+       else if(centrality<=45)  mycent=8;
+       else if(centrality<=50)  mycent=9;
+       else if(centrality<=55)  mycent=10;
+
+      else   mycent=11;
+
+  */
+
+ if(jentry%1000==0) cout << jentry<<endl;// event counter
+
+
+
+	float sumQxy[3][7][2] = {{{0}}}; //[ith][eta][x,y]
+	float multQv[7] = {0}; //[eta]
+
+
+ 
+hbimp->Fill(bimp);
+
+double refMult1 = 0;
+double refMult2 = 0;
+
+ 
+
+double Nch_L = 0;
+double Nch_R = 0;
+
+double Nch_L2 = 0;
+double Nch_R2 = 0;
+ 
+  for(int itrk=0;itrk<nh;itrk++) {  //track loop
+
+        float pt  = sqrt( TMath::Power(momx[itrk], 2.0 ) + TMath::Power(momy[itrk], 2.0 ) );
+	
+	float oldphi = atan2( momx[itrk], momy[itrk] );
+	float phi=oldphi;
+	
+	float the = atan2( pt, momz[itrk] );//atan2(pt/pz)
+	float eta = -log( tan( 0.5 * the ) );
+
+	//	phi += phiRP;
+	// float px = pt * cos(phi);
+	//	float py = pt * sin(phi);
+	//	pt  = sqrt( TMath::Power(px, 2.0 ) + TMath::Power(py, 2.0 ) );
+
+            h2pt->Fill(pt);
+            h2eta->Fill(eta);
+            h2phi->Fill(oldphi);
+            h2phis->Fill(phi);
+
+	    
+        if( pt>0.1 && fabs(eta)<0.5 ) refMult1++;
+	if( pt>0.0 && fabs(eta)<0.5 ) refMult2++;
+        if(eta >=  3.1 && eta <=  4.0) Nch_R++;
+        if(eta >= -4.0 && eta <= -3.1) Nch_L++;
+
+       
+        if( pt<0.15  || pt>2. ) continue;
+
+        int fEta = -1;
+
+	// TPC plane
+	if( eta>-1 && eta<-0.1 ) fEta = 0; // TPC East
+	if( eta>0.1 && eta<1 )   fEta = 1; // TPC West
+	if( fabs(eta)<0.1 )      fEta = 2; // TPC Mid
+
+      	// RXN plane
+       	if( eta>-3.0 && eta<-1.0 )   fEta = 3; // RXN East
+       	if( eta>1.0 && eta<3.0 )     fEta = 4; // RXN West
+
+	// BBC plane 
+	if( eta>-5 && eta<-3 )      fEta = 5; //East
+	if( eta>3.0 && eta<5  )     fEta = 6; //West
+
+	// if( fabs(eta)>1.1 && fabs(eta)<2.9 )     fEta = 7; // RXN combined
+        //if( fabs(eta)>3.0 && fabs(eta)<5.0 )     fEta = 8; // BBC combined
+	
+
+
+           if( fEta>-1 ){
+			for( int ith=0; ith<3; ith++ ){
+				sumQxy[ith][fEta][0] += pt * cos( (ith+2.0) * phi );
+				sumQxy[ith][fEta][1] += pt * sin( (ith+2.0) * phi );
+			}
+			multQv[fEta]++;
+	   }// end of eta selection
+      
+
+   }// end of track loop
+
+
+  
+      
+      if(Nch_L >= 1 && Nch_R >= 1){
+          href1 -> Fill(Nch_L + Nch_R);
+      }
+
+      if(Nch_L >= 2 && Nch_R >= 2){
+          href2 -> Fill(Nch_L + Nch_R);
+	  hbimp2->Fill(bimp);
+      }
+      if(refMult1>0){
+          href3 -> Fill(refMult1);
+	  hbimp3->Fill(bimp);
+      }
+      
+
+      if(refMult2>0){
+          href4 -> Fill(refMult1);
+	  hbimp4->Fill(bimp);
+      }
+
+      if(refMult1>2){
+          href5 -> Fill(refMult1);
+	  hbimp5->Fill(bimp);
+      }
+
+      float sumLR=Nch_L + Nch_R;
+  
+      // int fCent   = GetCentrality10_RefMult( refMult1 );// STAR def
+      //if( fCent<0 ) cout << fCent << endl;
+
+      // int fCent   =  GetCentrality10_RefMultPHENIX(sumLR);
+
+      int fCent   = GetCentrality10_Bimp(bimp);
+      //int fCent   =  GetCentrality10_BimpExp(bimp);
+      
+  	float fEP[3][7]; //[ith][eta]
+	float fQv[3][7]; //[ith][eta]
+	for( int ith=0; ith<3; ith++ ){ // flow harmonic loop
+	  for( int ieta=0; ieta<7; ieta++ ){ // ep detector gap 
+	    if( multQv[ieta]>5 ){ // multiplicity > 5
+				fEP[ith][ieta] = atan2( sumQxy[ith][ieta][1], sumQxy[ith][ieta][0] ) / ( ith + 2.0 );
+				fEP[ith][ieta] = atan2( sin( (ith+2.0)*fEP[ith][ieta] ), cos( (ith+2.0)*fEP[ith][ieta] ) );
+				fEP[ith][ieta] /= ( ith + 2.0 );
+
+ fQv[ith][ieta] = sqrt( TMath::Power( sumQxy[ith][ieta][0],2.0)+TMath::Power( sumQxy[ith][ieta][1], 2.0))/sqrt( multQv[ieta]);
+ H_Qw[ieta]->Fill( multQv[ieta] );
+			}else{
+				fEP[ith][ieta] = -9999;
+				fQv[ith][ieta] = -9999;
+	    }// end of mult cut selection
+	  } // end of loop on EP detectors
+	} // end of flow harmonic loop
+
+
+
+	for( int ith=0; ith<3; ith++ ){ // harmonic loop
+	  for( int ieta=0; ieta<7; ieta++ ){// eta EP detector loop
+	    if( fEP[ith][ieta]>-9000 ){ // EP reconstructed 
+				H_EP[ith][ieta]->Fill( fEP[ith][ieta] );
+				H_Qv[ith][ieta]->Fill( fQv[ith][ieta] );
+	    }// end of EP reconstructed
+	  }// end of eta loop
+	}// end of harm loop
+
+	
+		//Resolution
+	for( int ith=0; ith<3; ith++ ){
+		for( int icb=0; icb<3; icb++ ){
+			double psi1, psi2, fq1, fq2;
+       
+	       if ( icb==0 ){ psi1 = fEP[ith][0]; psi2 = fEP[ith][1]; fq1 = fQv[ith][0]; fq2 = fQv[ith][1]; } // TPC.E-TPC.W
+	       else if( icb==1 ){ psi1 = fEP[ith][3]; psi2 = fEP[ith][4]; fq1 = fQv[ith][3];  fq2 = fQv[ith][4]; } // RXN.E-RXN.W
+	       else { psi1 = fEP[ith][5]; psi2 = fEP[ith][6]; fq1 = fQv[ith][5]; fq2 = fQv[ith][6]; } // BBC.E-BBC.W
+		       
+
+		        
+		
+		
+
+			if( psi1<-9000 || psi2<-9000 ) continue;
+			if( fq1<0 || fq2<0 ) continue;
+
+
+
+                        double dPsi = ( ith + 2. ) * ( psi1 - psi2 );
+		   	dPsi = atan2( sin(dPsi), cos(dPsi) );
+			if(fCent>-1&&fCent<6){
+			
+			  HRes[ith][icb][fCent]->Fill(cos(dPsi) );
+
+			}
+
+		}
+	}
+
+
+
+// refmult star
+//float res2tpc[6]={0.503463,0.71591,0.749962,0.708934,0.61689,0.475386};
+//float res2rxn[6]={0.547883,0.791309,0.824213,0.793991,0.709638,0.561604};
+//float res2bbc[6]={0.252904,0.377894,0.401809,0.366931,0.300508,0.223784};
+
+// phenix ala cent
+
+// float res2tpc[6]={0.482228,0.704061,0.749667,0.7259,0.655706,0.541163};
+//float res2rxn[6]={ 0.526791,0.780304,0.824815,0.806745,0.746894,0.634004};
+//float res2bbc[6]={0.241286,0.375359,0.40532,0.384438,0.326188,0.250811};
+
+// 11.5 gev
+float res2tpc[6] = {0.323761,0.434068,0.433521,0.383028,0.307532,0.250556};
+float res2rxn[6] = {0.212259,0.333732,0.340171,0.303697,0.247778,0.206723};
+
+float res2bbc[6]={0.248559,0.383837,0.40206,0.373202,0.299496,0.2273};
+
+if(fCent>=0&&fCent<6){
+
+ 
+
+  
+  for(int itrk=0;itrk<nh;itrk++) {  //track loop
+
+        float pt  = sqrt( TMath::Power(momx[itrk], 2.0 ) + TMath::Power(momy[itrk], 2.0 ) );
+	
+	float oldphi = atan2( momx[itrk], momy[itrk] );
+	float phi=oldphi;
+	
+	float the = atan2( pt, momz[itrk] );//atan2(pt/pz)
+	float eta = -log( tan( 0.5 * the ) );
+
+        // phi += phiRP;
+	//  float px = pt * cos(phi);
+	//	float py = pt * sin(phi);
+	//	pt  = sqrt( TMath::Power(px, 2.0 ) + TMath::Power(py, 2.0 ) );
+
+
+	
+
+        if( pt<0.2  || pt>4.0 ) continue;
+
+
+	
+       int ipt = 0;
+ 	for(int i=0; i<npt-1;i++){
+	if(pt>=bin_w[i] && pt<bin_w[i+1]) ipt = i;
+	}// end of ipt loop
+
+        float v2tpc=-999.0;
+	float v2rxn=-999.0;
+	float v2bbc=-999.0;
+
+	if(eta>0&&eta<1.0){
+         v2tpc = cos(2.0 * (phi-fEP[0][0]) )/res2tpc[fCent];
+         v2rxn = cos(2.0 * (phi-fEP[0][3]) )/res2rxn[fCent];
+         v2bbc = cos(2.0 * (phi-fEP[0][5]) )/res2bbc[fCent];
+ 
+	}
+
+        if(eta<0&&eta>-1.0){
+         v2tpc = cos(2.0 * (phi-fEP[0][1]) )/res2tpc[fCent];
+         v2rxn = cos(2.0 * (phi-fEP[0][4]) )/res2rxn[fCent];
+         v2bbc = cos(2.0 * (phi-fEP[0][6]) )/res2bbc[fCent];
+ 
+	}
+
+        
+	
+
+        if(fabs(eta)<1.0){
+	  hv2[0][fCent][ipt][0]->Fill(v2tpc);
+          hv2[1][fCent][ipt][0]->Fill(v2rxn);
+          hv2[2][fCent][ipt][0]->Fill(v2bbc);
+	  hpt[fCent][ipt][0]->Fill(pt);
+
+          if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][0]->Fill(v2tpc);
+          hv22[1][ipt][0]->Fill(v2rxn);
+          hv22[2][ipt][0]->Fill(v2bbc);
+
+	  }
+	  
+
+	  if ( pdg[itrk]==211){
+              hv2[0][fCent][ipt][1]->Fill(v2tpc);
+              hv2[1][fCent][ipt][1]->Fill(v2rxn);
+              hv2[2][fCent][ipt][1]->Fill(v2bbc);
+	      hpt[fCent][ipt][1]->Fill(pt);
+
+
+          if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][1]->Fill(v2tpc);
+          hv22[1][ipt][1]->Fill(v2rxn);
+          hv22[2][ipt][1]->Fill(v2bbc);
+
+	  }
+
+
+	      
+
+          }// end of pion selection
+
+
+	  if ( pdg[itrk]==321){
+              hv2[0][fCent][ipt][2]->Fill(v2tpc);
+              hv2[1][fCent][ipt][2]->Fill(v2rxn);
+              hv2[2][fCent][ipt][2]->Fill(v2bbc);
+	      hpt[fCent][ipt][2]->Fill(pt);
+
+          if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][2]->Fill(v2tpc);
+          hv22[1][ipt][2]->Fill(v2rxn);
+          hv22[2][ipt][2]->Fill(v2bbc);
+
+	  }
+	      
+
+          }// end of kaon selection
+
+
+	   if ( pdg[itrk]==2212){
+              hv2[0][fCent][ipt][3]->Fill(v2tpc);
+              hv2[1][fCent][ipt][3]->Fill(v2rxn);
+              hv2[2][fCent][ipt][3]->Fill(v2bbc);
+	      hpt[fCent][ipt][3]->Fill(pt);
+
+           if(fCent>0&&fCent<4){
+
+          hv22[0][ipt][3]->Fill(v2tpc);
+          hv22[1][ipt][3]->Fill(v2rxn);
+          hv22[2][ipt][3]->Fill(v2bbc);
+
+	  }
+	      
+
+          }// end of proton selection
+
+	   
+	}
+
+        
+
+
+
+
+
+	
+
+
+  }// end of the track loop
+
+ }// end of centrality selection 
+
+
+
+
+
+
+
+
+
+	
+     
+ }
+
+
+
+// Urqmd 7.7 GeV
+
+int FlowANA::GetCentrality10_RefMult( double refMult ){
+
+	int fcent;
+	if     ( refMult>=268 ) fcent = 0; // 0-10% 
+	else if( refMult>=185 ) fcent = 1; //10-20% 
+	else if( refMult>=126 ) fcent = 2; //20-30% 
+	else if( refMult>=83 ) fcent = 3; //30-40% 
+	else if( refMult>=52 ) fcent = 4; //40-50% 
+	else if( refMult>= 21 ) fcent = 5; //50-60% 
+	else if( refMult>= 17 ) fcent = 6; //60-70% 
+	else if( refMult>= 9 ) fcent = 7; //70-80% 
+	else                    fcent =-1;
+
+	return fcent;
+}
+
+
+int FlowANA::GetCentrality10_RefMultPHENIX( double refMult ){
+
+	int fcent;
+	if     ( refMult>=516 ) fcent = 0; // 0-10% 
+	else if( refMult>=382 ) fcent = 1; //10-20% 
+	else if( refMult>=280 ) fcent = 2; //20-30% 
+	else if( refMult>=199 ) fcent = 3; //30-40% 
+	else if( refMult>=137 ) fcent = 4; //40-50% 
+	else if( refMult>= 90 ) fcent = 5; //50-60% 
+	else if( refMult>= 55 ) fcent = 6; //60-70% 
+	else if( refMult>= 31 ) fcent = 7; //70-80% 
+	else                    fcent =-1;
+
+	return fcent;
+}
+
+
+
+
+
+
+
+int FlowANA::GetCentrality10_Bimp( float bimp ){
+
+	int fcent;
+	if     ( bimp<4.5 ) fcent = 0; // 0-10%
+	else if( bimp<6.3 )  fcent = 1; //10-20%
+	else if( bimp<7.73 ) fcent = 2; //20-30%
+	else if( bimp<8.92 ) fcent = 3; //30-40%
+	else if( bimp<9.99)  fcent = 4; //40-50%
+	else if( bimp<10.83) fcent = 5; //50-60%
+	else if( bimp<11.78) fcent = 6; //60-70%
+	else if( bimp<12.61) fcent = 7; //70-80%
+	else                fcent =-1;
+
+	return fcent;
+}
+
+
+int FlowANA::GetCentrality10_BimpExp( float bimp ){
+
+	int fcent;
+	if     ( bimp<4.68 ) fcent = 0; // 0-10%
+	else if( bimp<6.47 ) fcent = 1; //10-20%
+	else if( bimp<7.99 ) fcent = 2; //20-30%
+	else if( bimp<9.31 ) fcent = 3; //30-40%
+	else if( bimp<10.48 ) fcent = 4; //40-50%
+	else if( bimp<11.49 ) fcent = 5; //50-60%
+	else if( bimp<12.35) fcent = 6; //60-70%
+	else if( bimp<13.0) fcent = 7; //70-80%
+	else                fcent =-1;
+
+	return fcent;
+}

Разница между файлами не показана из-за своего большого размера
+ 1168 - 0
FlowANA_test.h


+ 200 - 0
ReadFlow.C

@@ -0,0 +1,200 @@
+struct coord
+{
+    float xlow;
+    float ylow;
+    float xup;
+    float yup;
+};
+
+void ReadFlow(const char *inFileName)
+{
+    // Setting up global variables for the plot
+    gROOT->SetStyle("Pub");
+    gROOT->ForceStyle();
+    gStyle->SetPalette(kDarkRainBow);
+    gStyle->SetErrorX(0);
+    gStyle->SetTitleSize(0.05);
+
+    const double binpt [11] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.8,2.3,2.8,4.0};
+    const TString hName [4] = {TString("h^{#pm}"),TString("#pi^{#pm}"),TString("K^{#pm}"),TString("p+#bar{p}")};
+    const TString hCent [6] = {TString("0-10%"),TString("10-20%"),TString("20-30%"),TString("30-40%"),TString("40-50%"),TString("50-60%")};
+    const std::pair<Double_t,Double_t> v2range[6][4] = {{{-0.023,0.06},{-0.023,0.06},{-0.023,0.06},{-0.023,0.06}},
+                                                        {{-0.005,0.09},{-0.005,0.09},{-0.005,0.09},{-0.005,0.09}},
+                                                        {{-0.023,0.09},{-0.023,0.09},{-0.023,0.09},{-0.023,0.09}},
+                                                        {{-0.051,0.16},{-0.051,0.16},{-0.051,0.16},{-0.051,0.16}},
+                                                        {{-0.023,0.16},{-0.023,0.16},{-0.023,0.16},{-0.023,0.16}},
+                                                        {{-0.023,0.16},{-0.023,0.16},{-0.023,0.16},{-0.023,0.16}}};
+
+    TFile *fi = new TFile(inFileName,"read");
+    TH1F *hv2[4][6][10][4]; //idet-icent-ipt-ipid
+    TH1F *hv22[4][10][4];
+    TH1F *hpt[6][10][4];
+
+    for (int idet=0;idet<4;idet++)
+    {
+        for (int icent=0;icent<6;icent++)
+        {
+            for (int ipt=0; ipt<10; ipt++)
+            {
+                for (int ipid=0;ipid<4;ipid++)
+                {
+                    if (idet == 0)
+                    {
+                        hpt[icent][ipt][ipid] = (TH1F*) fi->Get(Form("hpt_%i_%i_%i",icent,ipt,ipid));
+                    }
+                    if (icent==0)
+                    {
+                        hv22[idet][ipt][ipid] = (TH1F*) fi->Get(Form("hv22_%i_%i_%i",idet,ipt,ipid));
+                    }
+                    hv2[idet][icent][ipt][ipid] = (TH1F*) fi->Get(Form("hv2_%i_%i_%i_%i",idet,icent,ipt,ipid));
+                }
+            }
+        }
+    }
+    std::vector<Double_t> vPt, vV2tpc, vV2fhcal;
+    std::vector<Double_t> ePt, eV2tpc, eV2fhcal;
+    int det1 = 0, det2 = 3;
+    TGraphErrors *grv2tpc[6][4], *grv2fhcal[6][4]; //icent-ipid
+
+    for (int icent=0;icent<6;icent++)
+    {
+        for (int ipid=0;ipid<4;ipid++)
+        {
+            for (int ipt=0;ipt<10;ipt++)
+            {
+                vPt.push_back(hpt[icent][ipt][ipid]->GetMean());
+                ePt.push_back(0.);
+                vV2tpc.push_back(hv2[det1][icent][ipt][ipid]->GetMean());
+                eV2tpc.push_back(hv2[det1][icent][ipt][ipid]->GetMeanError());
+                vV2fhcal.push_back(hv2[det2][icent][ipt][ipid]->GetMean());
+                eV2fhcal.push_back(hv2[det2][icent][ipt][ipid]->GetMeanError());
+            }
+            grv2tpc[icent][ipid] = new TGraphErrors(vPt.size(),&vPt[0],&vV2tpc[0],&ePt[0],&eV2tpc[0]);
+            grv2fhcal[icent][ipid] = new TGraphErrors(vPt.size(),&vPt[0],&vV2fhcal[0],&ePt[0],&eV2fhcal[0]);
+
+            grv2tpc[icent][ipid]->GetXaxis()->SetTitle("p_{T}, [GeV/c]");
+            grv2tpc[icent][ipid]->GetYaxis()->SetTitle("v_{2}");
+            grv2fhcal[icent][ipid]->GetXaxis()->SetTitle("p_{T}, [GeV/c]");
+            grv2fhcal[icent][ipid]->GetYaxis()->SetTitle("v_{2}");
+
+            grv2tpc[icent][ipid]  ->SetMarkerStyle(20);
+            grv2fhcal[icent][ipid]->SetMarkerStyle(26);
+            
+            grv2tpc[icent][ipid]  ->GetYaxis()->SetRangeUser(v2range[icent][ipid].first,v2range[icent][ipid].second);
+            grv2fhcal[icent][ipid]->GetYaxis()->SetRangeUser(v2range[icent][ipid].first,v2range[icent][ipid].second);
+
+            grv2tpc[icent][ipid]  ->GetXaxis()->SetLimits(0.,3.49);
+            grv2fhcal[icent][ipid]->GetXaxis()->SetLimits(0.,3.49);
+
+            grv2tpc[icent][ipid]  ->SetMarkerSize(1.6);
+            grv2fhcal[icent][ipid]->SetMarkerSize(1.6);
+
+            vPt.clear(); ePt.clear();
+            vV2tpc.clear(); eV2tpc.clear();
+            vV2fhcal.clear(); eV2fhcal.clear();
+        }
+    }
+
+    TH1F *hptTmp;
+    TGraphErrors *grv22tpc[4], *grv22fhcal[4];
+    for (int ipid=0;ipid<4;ipid++) {
+        for (int ipt = 0; ipt < 10; ipt++) {
+            hptTmp = (TH1F*) hpt[1][ipt][ipid]->Clone();
+            hptTmp->Add(hpt[2][ipt][ipid]);
+            hptTmp->Add(hpt[3][ipt][ipid]);
+            vPt.push_back(hptTmp->GetMean());
+            ePt.push_back(0.);
+            vV2tpc.push_back(hv22[det1][ipt][ipid]->GetMean());
+            eV2tpc.push_back(hv22[det1][ipt][ipid]->GetMeanError());
+            vV2fhcal.push_back(hv22[det2][ipt][ipid]->GetMean());
+            eV2fhcal.push_back(hv22[det2][ipt][ipid]->GetMeanError());
+            delete hptTmp;
+        }
+        grv22tpc[ipid] = new TGraphErrors(vPt.size(),&vPt[0],&vV2tpc[0],&ePt[0],&eV2tpc[0]);
+        grv22fhcal[ipid] = new TGraphErrors(vPt.size(),&vPt[0],&vV2fhcal[0],&ePt[0],&eV2fhcal[0]);
+
+        grv22tpc[ipid]->GetXaxis()->SetTitle("p_{T}, [GeV/c]");
+        grv22tpc[ipid]->GetYaxis()->SetTitle("v_{2}");
+        grv22fhcal[ipid]->GetXaxis()->SetTitle("p_{T}, [GeV/c]");
+        grv22fhcal[ipid]->GetYaxis()->SetTitle("v_{2}");
+
+        grv22tpc[ipid]  ->SetMarkerStyle(20);
+        grv22fhcal[ipid]->SetMarkerStyle(26);
+
+        grv22tpc[ipid]  ->GetYaxis()->SetRangeUser(v2range[1][ipid].first,v2range[1][ipid].second);
+        grv22fhcal[ipid]->GetYaxis()->SetRangeUser(v2range[1][ipid].first,v2range[1][ipid].second);
+
+        grv22tpc[ipid]  ->GetXaxis()->SetLimits(0.,3.49);
+        grv22fhcal[ipid]->GetXaxis()->SetLimits(0.,3.49);
+
+        grv22tpc[ipid]  ->SetMarkerSize(1.6);
+        grv22fhcal[ipid]->SetMarkerSize(1.6);
+
+        vPt.clear(); ePt.clear();
+        vV2tpc.clear(); eV2tpc.clear();
+        vV2fhcal.clear(); eV2fhcal.clear();
+    }
+
+    TCanvas *canv[6];
+    TLegend *leg[6][4];
+    TLine   *line0 = new TLine();
+    coord lc[4] = {{0.2,0.7,0.55,0.89},{0.2,0.82,0.55,0.89},{0.2,0.82,0.55,0.89},{0.2,0.82,0.55,0.89}};
+    line0->SetLineStyle(2);
+    line0->SetLineColor(4);
+    line0->SetLineWidth(1);
+    for (int icent=0; icent<6; icent++)
+    {
+        canv[icent] = new TCanvas(Form("canv_%i",icent),Form("canv cent %i",icent),900,900);
+        canv[icent] ->Divide(2,2);
+    }
+    for (int icent=0; icent<6; icent++)
+    {
+        for (int ipid=0;ipid<4;ipid++)
+        {
+            canv[icent]->cd(ipid+1);
+            grv2tpc[icent][ipid]->Draw("AP PLC PMC");
+            grv2fhcal[icent][ipid]->Draw("P PLC PMC");
+            leg[icent][ipid] = new TLegend(lc[ipid].xlow,lc[ipid].ylow,lc[ipid].xup,lc[ipid].yup);
+            leg[icent][ipid]->SetBorderSize(0.);
+            leg[icent][ipid]->SetHeader(Form("%s, %s",hCent[icent].Data(),hName[ipid].Data()),"C");
+            if (ipid == 0)
+            {
+                leg[icent][ipid]->AddEntry(grv2tpc[icent][ipid],"v2{TPC}","p");
+                leg[icent][ipid]->AddEntry(grv2fhcal[icent][ipid],"v2{FHCal}","p");
+            }
+            leg[icent][ipid]->Draw();
+
+            line0->DrawLine(0.,0.,3.49,0.);
+        }
+        // canv[icent]->SaveAs(Form("/home/peter/Documents/WorkLocal/STAR/Pics/Models/EPflow/v2_cent%i.png",icent));
+            if (icent == 0)
+                canv[icent]->Print(Form("/home/peter/Documents/WorkLocal/STAR/Pics/Models/EPflow/v2.pdf("),"pdf");
+            if (icent == 5)
+                canv[icent]->Print(Form("/home/peter/Documents/WorkLocal/STAR/Pics/Models/EPflow/v2.pdf)"),"pdf");
+            if (icent != 0 && icent != 5)
+                canv[icent]->Print(Form("/home/peter/Documents/WorkLocal/STAR/Pics/Models/EPflow/v2.pdf"),"pdf");
+    }
+
+    TCanvas *canv2 = new TCanvas("canv2","10-40%", 900, 900);
+    TLegend *leg2[4];
+    canv2->Divide(2,2);
+    for (int ipid=0;ipid<4;ipid++) {
+        canv2->cd(ipid + 1);
+        grv22tpc[ipid]->Draw("AP PLC PMC");
+        grv22fhcal[ipid]->Draw("P PLC PMC");
+
+        leg2[ipid] = new TLegend(lc[ipid].xlow,lc[ipid].ylow,lc[ipid].xup,lc[ipid].yup);
+        leg2[ipid]->SetBorderSize(0.);
+        leg2[ipid]->SetHeader(Form("10-40%, %s",hName[ipid].Data()),"C");
+        if (ipid == 0)
+        {
+            leg2[ipid]->AddEntry(grv22tpc[ipid],"v2{TPC}","p");
+            leg2[ipid]->AddEntry(grv22fhcal[ipid],"v2{FHCal}","p");
+        }
+        leg2[ipid]->Draw();
+
+        line0->DrawLine(0.,0.,3.49,0.);
+    }
+    canv2->SaveAs("/home/peter/Documents/WorkLocal/STAR/Pics/Models/EPflow/v2_cent1040.pdf");
+    
+}

+ 87 - 0
ReadRes.C

@@ -0,0 +1,87 @@
+#include <Math/SpecFuncMathMore.h>
+#include <TMath.h>
+
+R__LOAD_LIBRARY(libMathMore.so);
+
+// Resolution calculation
+//S----------------------------------------------------------------
+double GetRes(double _chi, double _harm){
+  double con = TMath::Sqrt(TMath::Pi() / 2) / 2;
+  double arg = _chi * _chi / 4.;
+  double order1 = (_harm - 1) / 2.;
+  double order2 = (_harm + 1) / 2.;
+  double res = con * _chi * exp(-arg) * (ROOT::Math::cyl_bessel_i(order1, arg) + ROOT::Math::cyl_bessel_i(order2, arg));
+  return res;
+}
+//E----------------------------------------------------------------
+
+
+// Chi2 calculation
+//S----------------------------------------------------------------
+double GetChi(double _res, double _harm, int accuracy){
+  double chi = 2.0;
+  double delta = 1.0;
+  for (int i = 0; i < accuracy; i++){
+    if(GetRes(chi, _harm) < _res) chi = chi + delta;
+    else chi = chi - delta;
+    delta = delta / 2.;
+  }
+  return chi;
+}
+//E----------------------------------------------------------------
+
+void ReadRes(const char *inFileName)
+{
+    TFile *fi = new TFile(inFileName,"read");
+    TH1F *hRes[3][4][6];
+    for (int ih=0; ih<3;ih++)
+    {
+        for (int icb=0;icb<4;icb++)
+        {
+            for (int ic=0; ic<6;ic++)
+            {
+                hRes[ih][icb][ic] = (TH1F*) fi->Get(Form("HRes_%i_%i_%i",ih,icb,ic));
+            }
+        }
+    }
+
+    std::vector<Double_t> vX, vY, vEX, vEY;
+    int cb = 3;
+    int h  = 0;
+    double chi, chiF, res2, res, resF;
+    std::cout << "Calculating resolution:" << std::endl;
+    for (int ic=0; ic<6; ic++)
+    {
+        res2 = hRes[h][cb][ic]->GetMean();
+        res = (res2>0) ? TMath::Sqrt(res2) : 0.;
+        chi = GetChi(res,2.,50);
+        chiF = TMath::Sqrt(2)*chi;
+        resF = GetRes(chiF,2.);
+        vX.push_back(10.*ic+5.);
+        vEX.push_back(0.);
+        if (res == 0){
+            vY.push_back(0.);
+            vEY.push_back(0.);
+        }
+        vY.push_back(resF);
+        vEY.push_back( hRes[h][cb][ic]->GetMeanError() * (resF/res) );
+
+        std::cout << "ic " << ic << ": " << resF << std::endl;
+    }
+    std::cout << std::endl;
+
+    std::cout << "res2fhcal = {";
+    for (int ic=0;ic<5;ic++)
+    {
+        std::cout << vY.at(ic) << ",";
+    }
+    std::cout << vY.at(5) << "};" << std::endl;
+
+    TGraphErrors *gr = new TGraphErrors(vX.size(),&vX[0],&vY[0],&vEX[0],&vEY[0]);
+    gr->SetMarkerStyle(21);
+    gr->SetMarkerSize(1.6);
+    gr->SetMarkerColor(1);
+    gr->SetLineColor(1);
+    
+    gr->Draw("AP");
+}

+ 28 - 0
main_proc.C

@@ -0,0 +1,28 @@
+//#include "FlowANA.C"
+#include "FlowANA.C"
+#include <TStopwatch.h>
+#include <TChain.h>
+#include <iostream>
+#include <fstream>
+
+void main_proc(const char *filelist, const char *output = "./test1.root")
+{
+  TStopwatch timer;
+  timer.Start();
+
+  TChain *chain = new TChain("mctree");
+  std::ifstream file(filelist);
+  std::string str;
+  while (std::getline(file,str))
+  {
+    chain->Add(str.c_str());
+  }
+
+  FlowANA *t = new FlowANA(chain);
+  t->ana_init(output);
+  t->Loop();
+  t->ana_end();
+
+  timer.Stop();
+  timer.Print();
+}

+ 28 - 0
main_proc_test.C

@@ -0,0 +1,28 @@
+//#include "FlowANA.C"
+#include "FlowANA_test.C"
+#include <TStopwatch.h>
+#include <TChain.h>
+#include <iostream>
+#include <fstream>
+
+void main_proc_test(const char *filelist, char *output)
+{
+  TStopwatch timer;
+  timer.Start();
+
+  TChain *chain = new TChain("mctree");
+  std::ifstream file(filelist);
+  std::string str;
+  while (std::getline(file,str))
+  {
+    chain->Add(str.c_str());
+  }
+
+  FlowANA *t = new FlowANA(chain);
+  t->ana_init(output);
+  t->Loop();
+  t->ana_end();
+
+  timer.Stop();
+  timer.Print();
+}

+ 252 - 0
res2.C

@@ -0,0 +1,252 @@
+#include "TFile.h"
+#include "TString.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+#include <iostream>
+#include "TLine.h"
+#include "TMath.h"
+#include <Math/SpecFuncMathMore.h>
+#include <TMath.h>
+
+#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
+R__LOAD_LIBRARY(libMathMore.so);
+#endif
+
+// Resolution calculation
+//S----------------------------------------------------------------
+double GetRes(double _chi, double _harm){
+  double con = TMath::Sqrt(TMath::Pi() / 2) / 2;
+  double arg = _chi * _chi / 4.;
+  double order1 = (_harm - 1) / 2.;
+  double order2 = (_harm + 1) / 2.;
+  double res = con * _chi * exp(-arg) * (ROOT::Math::cyl_bessel_i(order1, arg) + ROOT::Math::cyl_bessel_i(order2, arg));
+  return res;
+}
+//E----------------------------------------------------------------
+
+
+// Chi2 calculation
+//S----------------------------------------------------------------
+double GetChi(double _res, double _harm, int accuracy){
+  double chi = 2.0;
+  double delta = 1.0;
+  for (int i = 0; i < accuracy; i++){
+    if(GetRes(chi, _harm) < _res) chi = chi + delta;
+    else chi = chi - delta;
+    delta = delta / 2.;
+  }
+  return chi;
+}
+//E----------------------------------------------------------------
+
+
+
+
+void res2(Int_t mEnergy = 6)
+{
+
+#if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
+gSystem->Load(libMathMore.so);
+#endif
+
+//TFile *file = new TFile("sum.root");
+//TFile *file = new TFile("./out/res_11gev.root");
+TFile *file = new TFile("./out/res_7gev.root");
+  
+char hname1[800];
+char title[800];
+
+
+
+ 
+ float res2tpcEW[6];
+ float res2rxnEW[6]; 
+ float res2bbcEW[6];
+ float res2fhcalEW[6];
+
+float cent[6]={5,15,25,35,45,55};
+float eres[6];
+float ecent[6];
+
+float chi, res, chiF, resF;
+ 
+for (int ic=0; ic<6; ic++) {
+
+(void) sprintf(hname1,"HRes_%i_%i_%i",0,0,ic);
+TH1F *h3 = (TH1F*) file->Get(hname1);
+res = sqrt(h3->GetMean());
+chi = GetChi(res,1.,50);
+chiF = TMath::Sqrt(2)*chi;
+resF = GetRes(chiF,1.);
+res2tpcEW[ic]=resF;//sqrt(h3->GetMean());
+
+(void) sprintf(hname1,"HRes_%i_%i_%i",0,1,ic);
+TH1F *h4 = (TH1F*) file->Get(hname1);
+res2rxnEW[ic]=sqrt(h4->GetMean());
+
+(void) sprintf(hname1,"HRes_%i_%i_%i",0,2,ic);
+TH1F *h5 = (TH1F*) file->Get(hname1);
+res2bbcEW[ic]=sqrt(h5->GetMean());
+
+(void) sprintf(hname1,"HRes_%i_%i_%i",0,3,ic);
+TH1F *h6 = (TH1F*) file->Get(hname1);
+res = sqrt(h6->GetMean());
+chi = GetChi(res,2.,50);
+chiF = TMath::Sqrt(2)*chi;
+resF = GetRes(chiF,2.);
+res2fhcalEW[ic]=resF;
+
+ eres[ic]=0.01;
+ ecent[ic]=0.03;
+
+
+
+ //cout << res2bbcEW[ic] <<",";
+  
+ }
+ cout << endl;
+ cout << "Resolution for the analysis:" << endl;
+ cout << "float res2tpc[6] = {";
+ for (int ic=0; ic<6; ic++) {
+   cout << res2tpcEW[ic];
+   if (ic < 5) cout << ",";
+   else cout << "};" << endl;
+ }
+ cout << "float res2rxn[6] = {";
+ for (int ic=0; ic<6; ic++) {
+   cout << res2rxnEW[ic];
+   if (ic < 5) cout << ",";
+   else cout << "};" << endl;
+ }
+ cout << "float res2bbc[6] = {";
+ for (int ic=0; ic<6; ic++) {
+   cout << res2bbcEW[ic];
+   if (ic < 5) cout << ",";
+   else cout << "};" << endl;
+ }
+ cout << "float res2fhcal[6] = {";
+ for (int ic=0; ic<6; ic++) {
+   cout << res2fhcalEW[ic];
+   if (ic < 5) cout << ",";
+   else cout << "};" << endl;
+ }
+ cout << endl;
+
+ TCanvas *c1;
+ c1 = new TCanvas("c1","Flow analysis results",200,10,700,560);
+  
+ c1->GetFrame()->SetBorderSize(12);
+
+ c1->cd(); 
+	
+  gStyle->SetOptStat(0);
+
+ c1->GetFrame()->SetFillColor(21);
+ c1->GetFrame()->SetBorderSize(12);
+
+  gStyle->SetOptStat(0);
+  gStyle->SetPalette(0);
+  gStyle->SetCanvasColor(0);
+  gStyle->SetHistFillColor(10);
+  gStyle->SetHistFillStyle(0);
+  gStyle->SetOptTitle(1);
+  gStyle->SetOptStat(0);
+  c1->SetBorderMode(0);
+
+  c1->SetFillColor(0);
+  float xmin1=0.0;
+  float xmax1=60;
+  float ymin1=0.0;
+  float ymax1=1.0;
+
+
+
+sprintf(title,"PHSD: #sigma_{RP}   vs centrality for v_{2} , Au+Au  #sqrt{s_{NN}}=200 GeV");
+
+
+ TH2F *hr2 = new TH2F("hr2",title, 2,xmin1,xmax1,2,ymin1,ymax1);
+
+  hr2->SetXTitle("Centrality %");
+  hr2->SetYTitle("#sigma_{RP}");
+
+
+
+  hr2->Draw();
+
+
+  TGraphErrors *gr1;
+  TGraphErrors *gr2;
+  TGraphErrors *gr3;
+  TGraphErrors *gr4;
+
+
+
+
+  // const int npt=12;
+ gr1 = new TGraphErrors(6,cent,res2tpcEW,ecent,eres);
+ gr1->SetTitle("TGraphErrors Example ");
+ gr1->SetMarkerColor(kRed);
+ gr1->SetMarkerStyle(20);
+ gr1->SetMarkerSize(1.2);
+ gr1->Draw("P");
+
+
+
+  // const int npt=12;
+ gr2 = new TGraphErrors(6,cent,res2rxnEW,ecent,eres);
+ gr2->SetTitle("TGraphErrors Example ");
+ gr2->SetMarkerColor(kRed);
+ gr2->SetMarkerStyle(24);
+ gr2->SetMarkerSize(1.2);
+ gr2->Draw("P");
+
+
+
+
+  // const int npt=12;
+ gr3 = new TGraphErrors(6,cent,res2bbcEW,ecent,eres);
+ gr3->SetTitle("TGraphErrors Example ");
+ gr3->SetMarkerColor(kBlue);
+ gr3->SetMarkerStyle(21);
+ gr3->SetMarkerSize(1.2);
+ gr3->Draw("P");
+
+ // const int npt=12;
+ gr4 = new TGraphErrors(6,cent,res2fhcalEW,ecent,eres);
+ gr4->SetTitle("TGraphErrors Example ");
+ gr4->SetMarkerColor(kBlue);
+ gr4->SetMarkerStyle(25);
+ gr4->SetMarkerSize(1.2);
+ gr4->Draw("P");
+
+
+TF1 *f1 = new TF1("f1","[0]+[1]*(x)+[2]/x+[3]*x*x",5,55);
+gr3->Fit("f1","R+");
+gr3->Draw("P");
+
+
+
+ TLegend *legC12 = new TLegend(0.17,.75,0.27,.87);
+
+ legC12->SetFillColor(0);
+ legC12->SetBorderSize(0);
+ legC12->Draw();
+
+
+legC12->AddEntry(gr1,"TPC","lp");
+legC12->AddEntry(gr2,"RXN","lp");
+legC12->AddEntry(gr3,"BBC","lp");
+legC12->AddEntry(gr4,"FHCal","lp");
+
+
+
+
+ legC12->SetFillColor(0);
+ legC12->SetBorderSize(0);
+ legC12->Draw();
+
+
+
+
+
+}