void Point1_0() { TChain *ch1 = new TChain("tr"); for(Int_t i = 1; i<3; i++){ ch1->Add(TString::Format("AuAu200GeV_SM_NT3_DecaysOff_0_20fm_file0_%d.root", i)); } Int_t entries = ch1->GetEntries(); cout << entries << endl; Int_t Event, Mult; Float_t Imp; //Impact parameter ch1->SetBranchAddress("Event", &Event); ch1->SetBranchAddress("Mult", &Mult); ch1->SetBranchAddress("Imp", &Imp); Float_t Array[entries]; // Impact parameter array Int_t MArray[entries]; // Multiplicity array /* For multiplicities*/ Int_t PID[12000]; //Particle ID Float_t Mass[12000]; //Mass Float_t E[12000]; //Energy Float_t y[12000]; //Rapidity Float_t Pt[12000]; //Transverse Momentum Float_t P[120000]; // |P| mod P Float_t eta[12000]; //Pseudo-rapidity Float_t phi[12000]; //Azimuthal angle ch1->SetBranchAddress("PID",PID); ch1->SetBranchAddress("Mass",Mass); Float_t Px[12000]; Float_t Py[12000]; Float_t Pz[12000]; ch1->SetBranchAddress("Px",Px); ch1->SetBranchAddress("Py",Py); ch1->SetBranchAddress("Pz",Pz); /* For nucleons */ Int_t Stat; //Stat = 0 means Spectators Stat>0 means participants Float_t Nx; Float_t Ny; Float_t Nz; Float_t Nx_arr[394] = {}; Float_t Ny_arr[394] = {}; Float_t Nz_arr[394] = {}; Int_t Stat_arr[394] = {}; ch1->SetBranchAddress("Nx", &Nx_arr); ch1->SetBranchAddress("Ny", &Ny_arr); ch1->SetBranchAddress("Nz", &Nz_arr); ch1->SetBranchAddress("Stat", &Stat_arr); /* Rapidity (y), Pseudorapidity (eta), Transverse momentum (Pt), Azimuthal angle (phi) */ TCanvas *c8 = new TCanvas(); TH1F *hist1 = new TH1F("hist1", "Pt for S = 0", 100, 0, 4); for(Int_t i = 0; iGetEntry(i); for(Int_t j = 0; jFill(Pt[j]); } } c8->cd(); hist1->Draw(); /* Getting Collision Centrality from Impact Parameter Histogram */ TCanvas *c9 = new TCanvas(); TH1F *myhist = new TH1F("myhist", "Impact Parameter b", 100, 0., 0.); // auto binning for (Int_t i=0; iGetEntry(i); myhist->Fill(Imp); } c9->cd(); myhist->Draw(); /* Getting 0-5% centrality */ Int_t TBins = myhist->GetNbinsX(); // Total bins = 100 Int_t c = 0.05*entries; // Get 5% of Total events for 5% Centrality Int_t sum = 0; Int_t binnumber; Float_t bc = 0; // bin center value (Impact parameter) of a bin Double_t binentries[12000]; for(Int_t i = 0; iGetBinContent(i); sum = sum + binentries[i]; if(sum>c){ binnumber = i; bc = myhist->GetXaxis()->GetBinCenter(binnumber); cout << binnumber << " " << bc << " " << endl; break; } } /* Getting 70-80% centrality */ Int_t TBins_p = myhist->GetNbinsX(); Int_t p_1 = 0.7*entries; // 0-70% Int_t p_2 = 0.8*entries; // 0-80% Int_t sum_p; Int_t binnumber_p; Int_t B[2] = {p_1, p_2}; Float_t bc_p[2] = {}; // bin center value (Impact parameter) of a bin for 0-70% Double_t binentries_p[12000]; for(Int_t k = 0;k<2;k++){ sum_p = 0; for(Int_t l = 0; lGetBinContent(l); sum_p = sum_p + binentries_p[l]; if(sum_p>B[k]){ binnumber_p = l; bc_p[k] = myhist->GetXaxis()->GetBinCenter(binnumber_p); cout << binnumber_p << " " << bc_p[k] << " " << endl; break; } }} Float_t r_2[394] = {}; // r^2 = Nx^2 +Ny^2 Float_t r_2_sum[entries]; // r_2[0] + r_2[1] + .... + r_2[Mult]. It is for each collision. Index runs over i. r_2_sum[i] Float_t r_2_avg[entries]; // . It is for each collision. Index runs over i. r_2_avg[i] Float_t e2_cos[394] = {}; // cos term for each participant for n = 2 Float_t e2_sin[394] = {}; // sin term for each participant for n = 2 Float_t e2_cos_sum[entries]; // sum of cos terms for each participant for n = 2 Float_t e2_sin_sum[entries]; // sum of sin terms for each participant for n = 2 Float_t e2_cos_sum_avg[entries]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision Float_t e2_sin_sum_avg[entries]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision Float_t e2_num[entries]; // Numerator of eccentricity expression for n = 2 Float_t e2[entries]; //Eccentricity for each event for n = 2 Float_t e2_2[entries]; // Eccentricity^2 for each event for n = 2 Float_t e3_cos[394] = {}; // cos term for each participant for n = 3 Float_t e3_sin[394] = {}; // sin term for each participant for n = 3 Float_t e3_cos_sum[entries]; // sum of cos terms for each participant for n = 3 Float_t e3_sin_sum[entries]; // sum of sin terms for each participant for n = 3 Float_t e3_cos_sum_avg[entries]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision Float_t e3_sin_sum_avg[entries]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision Float_t e3_num[entries]; // Numerator of eccentricity expression for n = 3 Float_t e3[entries]; //Eccentricity for each event for n = 3 Float_t e3_2[entries]; // Eccentricity^2 for each event for n = 3 Float_t phi_part[394] = {}; //Azimuthal angle of participating nucleons Int_t N[entries]; // For counting the number of participants. Int_t cent = 0; Float_t E2[entries]; // Eccentricity array for n = 2 Float_t E3[entries]; // Eccentricity array for n = 3 Float_t diff_2_c[1] = {}; Float_t diff_3_c[1] = {}; Float_t err_2_c[1] = {}; Float_t err_3_c[1] = {}; Float_t error_2_c [1]= {}; Float_t error_3_c[1] = {}; Float_t sum_e2_2_c[1] = {}; Float_t sum_e3_2_c[1] = {}; Int_t counter_c[1] = {}; Float_t RMS_2_c[1] = {}; Float_t RMS_3_c[1] = {}; Float_t centr_c[1] = {2.5}; Float_t diff_2_p[1] = {}; Float_t diff_3_p[1] = {}; Float_t err_2_p[1] = {}; Float_t err_3_p[1] = {}; Float_t error_2_p[1]= {}; Float_t error_3_p[1] = {}; Float_t sum_e2_2_p[1] = {}; Float_t sum_e3_2_p[1] = {}; Int_t counter_p[1] = {}; Float_t RMS_2_p[1] = {}; Float_t RMS_3_p[1] = {}; Float_t centr_p[1] = {75}; for(Int_t i = 0; iGetEntry(i); Array[i] = Imp; MArray[i] = Mult; N[i] = 0; e2_2[i] = 0; e3_2[i] = 0; r_2_sum[i] = 0; e2_cos_sum[i] = 0; e2_sin_sum[i] = 0; e3_cos_sum[i] = 0; e3_sin_sum[i] = 0; for(Int_t j = 0; j<394; j++){ //cout << Nx_arr[i] << " " << Ny_arr[i] << " " << Nz_arr[i] << " " << Stat_arr[i] << endl; if(Stat_arr[j] > 0){ phi_part[j] = TMath::ATan2(Ny_arr[j],Nx_arr[j]); N[i]+=1; //counts the number of participating nucleons for each collision r_2[j] = ((TMath::Power(Nx_arr[j],2)) + (TMath::Power(Ny_arr[j],2))); r_2_sum[i] = r_2_sum[i] + r_2[j]; e2_cos[j]= r_2[j]*(TMath::Cos((2*phi_part[j]))); e2_sin[j]= r_2[j]*(TMath::Sin((2*phi_part[j]))); e2_cos_sum[i] = e2_cos_sum[i] + e2_cos[j]; e2_sin_sum[i] = e2_sin_sum[i] + e2_sin[j]; e3_cos[j]= r_2[j]*(TMath::Cos((3*phi_part[j]))); e3_sin[j]= r_2[j]*(TMath::Sin((3*phi_part[j]))); e3_cos_sum[i] = e3_cos_sum[i] + e3_cos[j]; e3_sin_sum[i] = e3_sin_sum[i] + e3_sin[j]; } } r_2_avg[i] = r_2_sum[i]/N[i]; e2_cos_sum_avg[i] = e2_cos_sum[i]/N[i]; e2_sin_sum_avg[i] = e2_sin_sum[i]/N[i]; e2_num[i] = TMath::Sqrt((TMath::Power(e2_cos_sum_avg[i],2)) + (TMath::Power(e2_sin_sum_avg[i],2))); //Numerator term for eccentricity e2[i] = e2_num[i]/r_2_avg[i]; e2_2[i] = TMath::Power(e2[i],2); e3_cos_sum_avg[i] = e3_cos_sum[i]/N[i]; e3_sin_sum_avg[i] = e3_sin_sum[i]/N[i]; e3_num[i] = TMath::Sqrt((TMath::Power(e3_cos_sum_avg[i],2)) + (TMath::Power(e3_sin_sum_avg[i],2))); //Numerator term for eccentricity e3[i] = e3_num[i]/r_2_avg[i]; e3_2[i] = TMath::Power(e3[i],2); E2[i] = e2[i]; E3[i] = e3[i]; if(Array[i]<=bc){ sum_e2_2_c[0] = sum_e2_2_c[0] + e2_2[i]; // 0 - 5% sum_e3_2_c[0] = sum_e3_2_c[0] + e3_2[i]; counter_c[0]+=1; } if(Array[i]>bc_p[0] && Array[i]<=bc_p[1]){ sum_e2_2_p[0] = sum_e2_2_p[0] + e2_2[i]; // 70 - 80% sum_e3_2_p[0] = sum_e3_2_p[0] + e3_2[i]; counter_p[0]+=1; } } RMS_2_c[0] = TMath::Sqrt((sum_e2_2_c[0]/counter_c[0])); RMS_3_c[0] = TMath::Sqrt((sum_e3_2_c[0]/counter_c[0])); RMS_2_p[0] = TMath::Sqrt((sum_e2_2_p[0]/counter_p[0])); RMS_3_p[0] = TMath::Sqrt((sum_e3_2_p[0]/counter_p[0])); cout << centr_c[0] << " " << RMS_2_c[0] << " " << RMS_3_c[0] << endl; cout << centr_p[0] << " " << RMS_2_p[0] << " " << RMS_3_p[0] << endl; for(Int_t i = 0; iGetEntry(i); //cout << E2[i] << endl; //cout << E3[i] << endl; if(Array[i]<=bc){ diff_2_c[0] = diff_2_c[0] + TMath::Power((RMS_2_c[0]-E2[i]),2); // 0 - 5% diff_3_c[0] = diff_3_c[0] + TMath::Power((RMS_3_c[0]-E3[i]),2); } else if(Array[i]>bc_p[0] && Array[i]<=bc_p[1]){ //cout << "bc_p[0] - bc_p[1]" << " " << RMS_2_p[1] << " " << E2[i] << endl; diff_2_p[0] = diff_2_p[0] + TMath::Power((RMS_2_p[0]-E2[i]),2); // 70 - 80% diff_3_p[0] = diff_3_p[0] + TMath::Power((RMS_3_p[0]-E3[i]),2); } } /* Calculating error for 0-5% centrality */ err_2_c[0] = TMath::Sqrt((diff_2_c[0]/(counter_c[0]-1))); err_3_c[0] = TMath::Sqrt((diff_3_c[0]/(counter_c[0]-1))); error_2_c[0] = err_2_c[0]/(TMath::Sqrt(counter_c[0])); error_3_c[0] = err_3_c[0]/(TMath::Sqrt(counter_c[0])); cout << error_2_c[0] << " " << error_3_c[0] << endl; /* Calculating error for 70-80% centrality */ err_2_p[0] = TMath::Sqrt((diff_2_p[0]/(counter_p[0]-1))); err_3_p[0] = TMath::Sqrt((diff_3_p[0]/(counter_p[0]-1))); error_2_p[0] = err_2_p[0]/(TMath::Sqrt(counter_p[0])); error_3_p[0] = err_3_p[0]/(TMath::Sqrt(counter_p[0])); cout << error_2_p[0] << " " << error_3_p[0] << endl; /* Graph with error 0-5% Centrality */ TGraphErrors* er1 = new TGraphErrors(); er1->SetPoint(1, centr_c[0], RMS_2_c[0]); er1->SetPointError(1, 0., error_2_c[0]); er1->SetTitle("RMS Eccentricity vs Centrality % (Imp b) n = 2 S = 0;Centrality %(Imp);RMS Eccentricity n = 2"); er1->SetLineColor(40); er1->SetMarkerColor(1); er1->SetMarkerStyle(20); TCanvas *c20 = new TCanvas(); c20->cd(); er1->Draw("AP"); TGraphErrors* er2 = new TGraphErrors(); er2->SetPoint(1, centr_c[0], RMS_3_c[0]); er2->SetPointError(1, 0., error_3_c[0]); er2->SetTitle("RMS Eccentricity vs Centrality % (Imp b) n = 3 S = 0;Centrality %(Imp);RMS Eccentricity n = 3"); er2->SetLineColor(41); er2->SetMarkerColor(2); er2->SetMarkerStyle(21); TCanvas *c21 = new TCanvas(); c21->cd(); er2->Draw("AP"); /* Graph with error 70-80% Centrality */ TGraphErrors* er3 = new TGraphErrors(); er3->SetPoint(1, centr_p[0], RMS_2_p[0]); er3->SetPointError(1, 0., error_2_p[0]); er3->SetTitle("RMS Eccentricity vs Centrality % (Imp b) n = 2 S = 0;Centrality %(Imp);RMS Eccentricity n = 2"); er3->SetLineColor(40); er3->SetMarkerColor(1); er3->SetMarkerStyle(20); TCanvas *c22 = new TCanvas(); c22->cd(); er3->Draw("AP"); TGraphErrors* er4 = new TGraphErrors(); er4->SetPoint(1, centr_p[0], RMS_3_p[0]); er4->SetPointError(1, 0., error_3_p[0]); er4->SetTitle("RMS Eccentricity vs Centrality % (Imp b) n = 3 S = 0;Centrality %(Imp);RMS Eccentricity n = 3"); er4->SetLineColor(41); er4->SetMarkerColor(2); er4->SetMarkerStyle(21); TCanvas *c23 = new TCanvas(); c23->cd(); er2->Draw("AP"); /*TCanvas *c22 = new TCanvas(); auto MG1 = new TMultiGraph(); MG1->SetTitle("RMS Eccentricity vs Centrality % (Imp b) for S = 0;Centrality % (Imp); RMS Eccentricity"); MG1->Add(er1, "cp"); MG1->Add(er2, "cp"); MG1->Draw("acp"); c22->BuildLegend(); TFile* f1 = new TFile("Central_0.root", "CREATE"); er1->Write("Graph1"); er2->Write("Graph2"); MG1->Write("GS0"); hist1->Write(); */ //input->Close(); }