Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:13:28 

0001 #include <TNtuple.h>
0002 #include <TF1.h>
0003 #include <TLine.h>
0004 #include "sPHENIXStyle/sPhenixStyle.C"
0005 
0006 #include <iostream>
0007 #include <fstream>
0008 
0009 //----- MY STUFF -------------------------------------------------
0010 double CBFunction(double *x, double *p)
0011 {
0012   double norm = p[0];
0013   double alpha = p[1];  // tail start
0014   double n = p[2];      // tail shape 
0015   double mu = p[3];     // upsilon mass
0016   double sigma = p[4];  // upsilon width
0017 
0018   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0019   double B = n/fabs(alpha) - fabs(alpha);
0020   double k = (x[0]-mu)/sigma;
0021 
0022   double val;
0023   if( k > -alpha )
0024     val = norm*TMath::Exp(-0.5*pow(k,2));
0025   else
0026     val = norm*A*pow(B-k,-n);
0027 
0028   if( TMath::IsNaN(val) )
0029     val = 0.0;
0030 
0031   return val;
0032 }
0033 
0034 double CBFunction_withBG(double *x, double *p)
0035 {
0036   double norm = p[0];
0037   double alpha = p[1];  // tail start
0038   double n = p[2];      // tail shape 
0039   double mu = p[3];     // upsilon mass
0040   double sigma = p[4];  // upsilon width
0041 
0042   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0043   double B = n/fabs(alpha) - fabs(alpha);
0044   double k = (x[0]-mu)/sigma;
0045 
0046   double val;
0047   if( k > -alpha )
0048     val = norm*TMath::Exp(-0.5*pow(k,2));
0049   else
0050     val = norm*A*pow(B-k,-n);
0051 
0052   double bgnorm1 = p[5];
0053   double bgslope1 = p[6];
0054   double bg = exp(bgnorm1+x[0]*bgslope1);
0055 
0056   val = val + bg;
0057 
0058   if( TMath::IsNaN(val) )
0059     val = 0.0;
0060 
0061   return val;
0062 }
0063 //---- ALICE STUFF ------------------------------------------------------------------
0064 
0065 Double_t CrystallBall2( Double_t x, Double_t mean, Double_t sigma, Double_t alpha1, Double_t n1, Double_t alpha2, Double_t n2 )
0066 {
0067   Double_t t = (x-mean)/sigma;
0068   if( t < -alpha1 )
0069   {
0070     Double_t a = TMath::Power( n1/alpha1, n1 )*TMath::Exp( -TMath::Power( alpha1, 2 )/2 );
0071     Double_t b = n1/alpha1 - alpha1;
0072     return a/TMath::Power( b - t, n1 );
0073   } else if( t > alpha2 ) {
0074 
0075     Double_t a = TMath::Power( n2/alpha2, n2 )*TMath::Exp( -TMath::Power( alpha2,2 )/2 );
0076     Double_t b = n2/alpha2 - alpha2;
0077     return a/TMath::Power( b + t, n2 );
0078   } else return TMath::Exp( -TMath::Power( t, 2 )/2 );
0079 }
0080 
0081 Double_t CrystallBall2Integral( Double_t sigma, Double_t alpha1, Double_t n1, Double_t alpha2, Double_t n2 )
0082 {
0083 // get corresponding integral
0084   alpha1 = fabs( alpha1 );
0085   alpha2 = fabs( alpha2 );
0086   return sigma*(
0087                 n1/(alpha1*(n1-1))*TMath::Exp( -pow( alpha1, 2 )/2 ) +
0088                 n2/(alpha2*(n2-1))*TMath::Exp( -pow( alpha2,2 )/2 ) +
0089     TMath::Sqrt( TMath::Pi()/2)*TMath::Erfc( -alpha1/TMath::Sqrt(2) ) -
0090     TMath::Sqrt( TMath::Pi()/2)*TMath::Erfc( alpha2/TMath::Sqrt(2) ) );
0091 }
0092 
0093 Double_t CrystallBall2( Double_t *x, Double_t *par )
0094 {
0095   // get normalized Crystal ball
0096   Double_t result = CrystallBall2( x[0], par[1], par[2], par[3], par[4], par[5], par[6] );
0097   // get integral
0098   Double_t integral = CrystallBall2Integral( par[2], par[3], par[4], par[5], par[6] );
0099   // return scaled Crystalball so that par[0] corresponds to integral
0100   return par[0] * result/integral;
0101 }
0102 
0103 //======================================================================
0104 
0105 void plotembed_vscent()
0106 {
0107 SetsPhenixStyle();
0108 gStyle->SetOptStat(0);
0109 gStyle->SetOptFit(1);
0110 
0111 double RR = 15.;
0112 const int nbins = 10;
0113 double bins[nbins+1];
0114 double AA = RR*RR/double(nbins);
0115 bins[0] = 0.;
0116 bins[nbins] = 16.;
0117 for(int i=1; i<nbins; i++) { bins[i] = sqrt(i*AA); }
0118 for(int i=0; i<=nbins; i++) { cout << "bin edge: " << i << " " << bins[i] << endl; }
0119 
0120 TF1* fCB = new TF1("fCB",CBFunction,6.,12.,5);
0121 TF1* fCBbg = new TF1("fCBbg",CBFunction_withBG,6.,12.,7);
0122 
0123 float mass,type,pt,eta,rap,pt1,pt2,eta1,eta2,bimp;
0124 float chisq1,chisq2,dca2d1,dca2d2,eop3x3_1,eop3x3_2;
0125 float nmvtx1,nmvtx2,ntpc1,ntpc2;
0126 
0127 char hname[99];
0128 TH1D* hm[nbins];
0129 TH1D* hmss[nbins];
0130 TH1D* hmnbg[nbins];
0131 TH1D* hbimp[nbins];
0132 for(int i=0; i<nbins; i++) {
0133   sprintf(hname,"hm_%d",i);
0134   hm[i] = new TH1D(hname,hname,120,6.,12.);
0135   sprintf(hname,"hmss_%d",i);
0136   hmss[i] = new TH1D(hname,hname,120,6.,12.);
0137   sprintf(hname,"hmnbg_%d",i);
0138   hmnbg[i] = new TH1D(hname,hname,120,6.,12.);
0139   sprintf(hname,"hbimp_%d",i);
0140   hbimp[i] = new TH1D(hname,hname,1000,0.,20.);
0141   hm[i]->Sumw2();
0142   hmss[i]->Sumw2();
0143   hmnbg[i]->Sumw2();
0144 }
0145 
0146 TH1D* hmass = new TH1D("hmass","",120,6.,12.);
0147 TH1D* hmass_ss = new TH1D("hmass_ss","",120,6.,12.);
0148 TH1D* hmass_nobg = new TH1D("hmass_nobg","",120,6.,12.);
0149 hmass->Sumw2();
0150 hmass_ss->Sumw2();
0151 hmass_nobg->Sumw2();
0152 
0153 TLine* l1 = new TLine(7.,0.,11.,0.);
0154 l1->SetLineStyle(2);
0155 
0156 //TNtuple* ntp2;
0157 TChain* cntp2 =  new TChain("ntp2");
0158 TChain* cntp2c =  new TChain("ntp2");
0159 
0160 //-------------------------------------------------------------------
0161 
0162 // central ---------------------------------------
0163 
0164 string infname_central;
0165 ifstream ifs_central;
0166 ifs_central.open("centrallist.txt");
0167 if ( ifs_central.fail() ) { cout << "FAIL TO READ INPUT FILE 1" << endl; ifs_central.close(); return; }
0168 while(!ifs_central.eof()) {
0169   ifs_central >> infname_central;
0170   cntp2c->AddFile(infname_central.c_str());
0171 }
0172 ifs_central.close();
0173 
0174 cout << "number of CENTRAL entries = " << cntp2c->GetEntries() << endl;
0175 
0176 cntp2c->SetBranchAddress("type",&type);
0177 cntp2c->SetBranchAddress("mass",&mass);
0178 cntp2c->SetBranchAddress("pt",&pt);
0179 cntp2c->SetBranchAddress("eta",&eta);
0180 cntp2c->SetBranchAddress("rap",&rap);
0181 cntp2c->SetBranchAddress("pt1",&pt1);
0182 cntp2c->SetBranchAddress("pt2",&pt2);
0183 cntp2c->SetBranchAddress("eta1",&eta1);
0184 cntp2c->SetBranchAddress("eta2",&eta2);
0185 cntp2c->SetBranchAddress("chisq1",&chisq1);
0186 cntp2c->SetBranchAddress("dca2d1",&dca2d1);
0187 cntp2c->SetBranchAddress("eop3x3_1",&eop3x3_1);
0188 cntp2c->SetBranchAddress("chisq2",&chisq2);
0189 cntp2c->SetBranchAddress("dca2d2",&dca2d2);
0190 cntp2c->SetBranchAddress("eop3x3_2",&eop3x3_2);
0191 cntp2c->SetBranchAddress("nmvtx1",&nmvtx1);
0192 cntp2c->SetBranchAddress("nmvtx2",&nmvtx2);
0193 cntp2c->SetBranchAddress("ntpc1",&ntpc1);
0194 cntp2c->SetBranchAddress("ntpc2",&ntpc2);
0195 
0196   for(int j=0; j<cntp2c->GetEntries(); j++) {
0197     cntp2c->GetEvent(j);
0198     if(chisq1>3 || chisq2>3) continue;
0199     if(nmvtx1<2 || nmvtx2<2) continue;
0200     if(ntpc1<26 || ntpc2<26) continue;
0201     if(eop3x3_1<0.7 || eop3x3_2<0.7) continue;
0202     //if(type==1) { hmass->Fill(mass); hmass_nobg->Fill(mass); }
0203     //if(type==2 || type==3) { hmass_ss->Fill(mass); }
0204       if(type==1) { hm[0]->Fill(mass); hmnbg[0]->Fill(mass); }
0205       if(type==2 || type==3) { hmss[0]->Fill(mass); }
0206   }
0207 // end central -----------------------------------
0208 
0209 // minbias ---------------------------------------
0210 string infname;
0211 ifstream ifs;
0212 ifs.open("mblist.txt");
0213 if ( ifs.fail() ) { cout << "FAIL TO READ INPUT FILE 1" << endl; ifs.close(); return; }
0214 while(!ifs.eof()) {
0215   ifs >> infname;
0216   cntp2->AddFile(infname.c_str());
0217 }
0218 ifs.close();
0219 
0220 cout << "number of MB entries = " << cntp2->GetEntries() << endl;
0221 
0222 cntp2->SetBranchAddress("b",&bimp);
0223 cntp2->SetBranchAddress("type",&type);
0224 cntp2->SetBranchAddress("mass",&mass);
0225 cntp2->SetBranchAddress("pt",&pt);
0226 cntp2->SetBranchAddress("eta",&eta);
0227 cntp2->SetBranchAddress("rap",&rap);
0228 cntp2->SetBranchAddress("pt1",&pt1);
0229 cntp2->SetBranchAddress("pt2",&pt2);
0230 cntp2->SetBranchAddress("eta1",&eta1);
0231 cntp2->SetBranchAddress("eta2",&eta2);
0232 cntp2->SetBranchAddress("chisq1",&chisq1);
0233 cntp2->SetBranchAddress("dca2d1",&dca2d1);
0234 cntp2->SetBranchAddress("eop3x3_1",&eop3x3_1);
0235 cntp2->SetBranchAddress("chisq2",&chisq2);
0236 cntp2->SetBranchAddress("dca2d2",&dca2d2);
0237 cntp2->SetBranchAddress("eop3x3_2",&eop3x3_2);
0238 cntp2->SetBranchAddress("nmvtx1",&nmvtx1);
0239 cntp2->SetBranchAddress("nmvtx2",&nmvtx2);
0240 cntp2->SetBranchAddress("ntpc1",&ntpc1);
0241 cntp2->SetBranchAddress("ntpc2",&ntpc2);
0242 
0243   for(int j=0; j<cntp2->GetEntries(); j++) {
0244     cntp2->GetEvent(j);
0245     if(chisq1>3 || chisq2>3) continue;
0246     if(nmvtx1<2 || nmvtx2<2) continue;
0247     if(ntpc1<26 || ntpc2<26) continue;
0248     if(eop3x3_1<0.7 || eop3x3_2<0.7) continue;
0249     if(type==1) { hmass->Fill(mass); hmass_nobg->Fill(mass); }
0250     if(type==2 || type==3) { hmass_ss->Fill(mass); }
0251 
0252     for(int j=0; j<nbins; j++) {
0253       if(bimp>bins[j] && bimp<bins[j+1]) {
0254         if(type==1) { hm[j]->Fill(mass); hmnbg[j]->Fill(mass); hbimp[j]->Fill(bimp); }
0255         if(type==2 || type==3) { hmss[j]->Fill(mass); }
0256       }
0257     }    
0258 
0259   }
0260 // end minbias -------------------------------------------------------
0261 
0262 for(int j=0; j<nbins; j++) { hmnbg[j]->Add(hmss[j],-1.0); }
0263 hmass_nobg->Add(hmass_ss,-1.0);
0264 
0265 double massres[nbins];
0266 double massreserr[nbins];
0267 double x[nbins];
0268 for(int i=0; i<nbins; i++) { x[i] = hbimp[i]->GetMean(); cout << "mean b = " << x[i] << endl; }
0269 
0270 //----------------------------------------------------------------
0271 
0272 // ALICE Double Crystal Ball function
0273 TF1 *f2 = new TF1("f2",CrystallBall2,7,11,7);
0274   f2->SetParameter(0, 2000. );
0275   f2->SetParameter(1, 9.46 );
0276   f2->SetParameter(2, 0.13 );
0277   f2->SetParameter(3, 1);
0278   f2->SetParameter(4, 3);
0279   f2->SetParameter(5, 1 );
0280   f2->SetParameter(6, 5 );
0281   f2->SetParNames("normalization", "mean", "sigma","alpha1","n1","alpha2","n2");
0282     f2->SetParLimits(1, 9.40, 9.55);
0283     f2->SetParLimits(2, 0.06, 0.20);
0284     f2->SetParLimits(3, 0.120, 10);
0285     f2->SetParLimits(4, 1.01, 10);
0286     f2->SetParLimits(5, 0.1, 10);
0287     f2->SetParLimits(6, 1.01, 10);
0288 
0289 
0290 TCanvas* c3 = new TCanvas("c3","inv. mass allpT",20,20,600,500);
0291 
0292 hmass->Draw();
0293 hmass_ss->SetLineColor(kRed);
0294 hmass_ss->Scale(1.2);
0295 hmass_ss->Draw("same");
0296 
0297 TCanvas* c33 = new TCanvas("c33","inv. mass nobg",20,20,600,500);
0298 //c33->Divide(3,2);
0299 /*
0300   fCB->SetParameter(0,1000.);
0301   fCB->SetParameter(1,1.0);
0302   fCB->SetParameter(2,1.0);
0303   fCB->SetParLimits(2,0.1,99.);
0304   fCB->SetParameter(3,9.0);
0305   fCB->SetParameter(4,0.100);
0306 */
0307 
0308 //  c33->cd(1);
0309   hmass_nobg->SetAxisRange(7.,11.);
0310   hmass_nobg->Fit("f2","rlq","",7.5,10.0);
0311 
0312   for(int j=0; j<nbins; j++) { 
0313 //    c33->cd(j+2);
0314   fCB->SetParameter(0,hmnbg[j]->GetMaximum());
0315     hmnbg[j]->SetTitle(";Invariant mass, GeV ; Counts");
0316     hmnbg[j]->Fit("f2","rlq","",7.5,10.0); 
0317     hmnbg[j]->SetAxisRange(7.,11.);
0318 //    l1->Draw();
0319     cout << "mass resolution = " << j << " " << f2->GetParameter(2) << " +- " << f2->GetParError(2) << endl;
0320     massres[j] = 1000.*f2->GetParameter(2);
0321     massreserr[j] = 1000.*f2->GetParError(2);
0322   }
0323   hmass_nobg->SetTitle(";Invariant mass, GeV ; Counts");
0324   hmass_nobg->Draw();
0325   TLatex* lat1 = new TLatex(7.2,hmass_nobg->GetMaximum()*0.4,"Minimum Bias"); lat1->Draw();
0326   l1->Draw();
0327 
0328 
0329 //  double fitndf = hmass_nobg->GetFunction("fCB")->GetNDF();
0330 //  double fitchi2 = hmass_nobg->GetFunction("fCB")->GetChisquare();
0331 //  cout << "chi2 = " << fitchi2/fitndf << endl;
0332 
0333 double scale = 1.4;
0334 TF1* fws = new TF1("fws","[0]/(1+exp((x-[1])/[2]))",0.,20.);
0335 fws->SetParameter(0,8.68398e+01*scale);
0336 fws->SetParameter(1,1.51845e+01);
0337 fws->SetParameter(2,6.52512e-01);
0338 TF1* fws_sum = new TF1("fws_sum","(x*[0])/(1+exp((x-[1])/[2]))",0.,20.);
0339 fws_sum->SetParameter(0,8.68398e+01*scale);
0340 fws_sum->SetParameter(1,1.51845e+01);
0341 fws_sum->SetParameter(2,6.52512e-01);
0342 double avg[nbins];
0343 for(int j=0; j<nbins; j++) {
0344   avg[j] = fws_sum->Integral(bins[j],bins[j+1]);
0345 //  cout << "expected = " << avg[j] << endl;
0346 }
0347 avg[0] = avg[0]*0.9;
0348 
0349 double sumstart = 8.8;
0350 double sumstop  = 9.8;
0351 int fbin = hmass_nobg->FindBin(sumstart + 0.001);
0352 int lbin = hmass_nobg->FindBin(sumstop  - 0.001);
0353 cout << "counting bins from " << fbin << " to " << lbin << endl;
0354 cout << "   in mass range from " << hmass_nobg->GetBinLowEdge(fbin) << " to " << hmass_nobg->GetBinLowEdge(lbin)+hmass_nobg->GetBinWidth(lbin) << endl;
0355 
0356   double sum =0.;
0357   double sumsum[nbins]; for(int i=0; i<nbins; i++) { sumsum[i]=0.; }
0358   double y[nbins],ey[nbins];
0359   for(int k=fbin; k<=lbin; k++) {
0360     sum += hmass_nobg->GetBinContent(k);
0361     for(int i=0; i<nbins; i++) { 
0362       sumsum[i] += hmnbg[i]->GetBinContent(k); 
0363       y[i] = sumsum[i]/avg[i];
0364       ey[i] = y[i]*0.05;
0365     }
0366   }
0367   cout << "Number of Upsilons = " << sum << endl;
0368   for(int i=0; i<nbins; i++) { cout << sumsum[i] << endl; }
0369 
0370   cout << "Number of Upsilons from fit = " << fCB->Integral(sumstart,sumstop)/hmass_nobg->GetBinWidth(lbin) << " " << fCB->Integral(5.0,11.0)/hmass_nobg->GetBinWidth(lbin) << endl;
0371 
0372 TCanvas* c4 = new TCanvas("c4","mass resolution vs Bimp",20,20,600,500);
0373 TH2D* hh = new TH2D("hh","",10,0.,18.,10,60.,140.);
0374 hh->SetTitle(";Impact parameter, fm ; Mass resolution, MeV");
0375 hh->Draw();
0376 
0377 TGraphErrors* gr1 = new TGraphErrors(nbins,x,massres,0,massreserr);
0378 gr1->SetMarkerStyle(20);
0379 gr1->SetMarkerColor(kBlack);
0380 gr1->SetMarkerSize(1.5);
0381 gr1->SetLineColor(kBlack);
0382 gr1->SetLineWidth(2);
0383 gr1->Draw("p");
0384 
0385 /*
0386 TCanvas* c44 = new TCanvas("c44","mass resolution vs Bimp",120,120,600,500);
0387 TH2D* hhh = new TH2D("hhh","",10,0.,18.,10,0.,1.5);
0388 hhh->Draw();
0389 
0390 TGraphErrors* gr2 = new TGraphErrors(nbins,x,y,0,ey);
0391 gr2->SetMarkerStyle(20);
0392 gr2->SetMarkerColor(kBlue);
0393 gr2->SetMarkerSize(1.5);
0394 gr2->SetLineColor(kBlue);
0395 gr2->Draw("p");
0396 */
0397 
0398 TCanvas* c5 = new TCanvas("c5","",60,60,600,500);
0399 hmnbg[0]->Draw();
0400 TLatex* lat2 = new TLatex(7.2,hmnbg[0]->GetMaximum()*0.4,"10% central Au+Au"); lat2->Draw();
0401 l1->Draw();
0402 
0403 }
0404
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!