本文整理汇总了C++中TH2D::Fill方法的典型用法代码示例。如果您正苦于以下问题：C++ TH2D::Fill方法的具体用法？C++ TH2D::Fill怎么用？C++ TH2D::Fill使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TH2D的用法示例。
在下文中一共展示了TH2D::Fill方法的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: mbfit
//.........这里部分代码省略.........
// wgt = wgt*2.0/tree1->GetLeaf("MBScalarEt")->GetValue(0);
//}
// / (0.5 * (1. + TMath::Erf((pmcsmbz-.5)/(TMath::Sqrt(2)*1.0))));
// + (0.5 * (1. + TMath::Erf((pmcsmbz-.5)/(TMath::Sqrt(2)*5)))); //0.5 is good for upper region, 0.75 for middle
//cout<<TMath::Power(pmcsmbz, 0.5/(10/(TMath::Sqrt(pmcsmbz))))<<" "<<(0.01 * (1. + TMath::Erf((pmcsmbz-5)/(TMath::Sqrt(2)))))<<endl;
//wgt = 0.582 * (1. + TMath::Erf((pmcsmbz-1.0935)/(TMath::Sqrt(2)*4.3111)));
// wgt = wgt*(0.53 * (1. + TMath::Erf((pmcsmbz-0.288)/(TMath::Sqrt(2)*4.0453))));
// wgt = wgt*(0.489 * (1. + TMath::Erf((pmcsmbz+0.305)/(TMath::Sqrt(2)*4.89144))));
/**
//works ok-ish but doesn't work as weight (obviously...)
pmcsmbz = tree1->GetLeaf("MBScalarEt")->GetValue(0)*1.5;
pmcsmbz = pmcsmbz*(0.544657 * (1. + TMath::Erf((pmcsmbz-0.6808)/(TMath::Sqrt(2)*3.18052))));
pmcsmbz = pmcsmbz*(0.530694 * (1. + TMath::Erf((pmcsmbz-0.726896)/(TMath::Sqrt(2)*1.0429))));
pmcsmbz = pmcsmbz*(0.480563 * (1. + TMath::Erf((pmcsmbz-0.321864)/(TMath::Sqrt(2)*0.536511))));
**/
// pmcsmbz = pmcsmbz*(0.342784 * (1. + TMath::Erf((pmcsmbz-0.167683)/(TMath::Sqrt(2)*0.540161))));
//pmcsmbz = pmcsmbz*(0.9587 + 0.00365*pmcsmbz);
// pmcsmbz = pmcsmbz*(0.903 + 0.007296*pmcsmbz + 0.00003*pmcsmbz*pmcsmbz - 0.00000119*pmcsmbz*pmcsmbz*pmcsmbz);//*ratiotest2->FindBin(tree1->GetLeaf("MBScalarEt")->GetValue(0));
// if (tree1->GetLeaf("ScalarEt")->GetValue(0) > 30) pmcsmbz = tree1->GetLeaf("MBScalarEt")->GetValue(0)*1.0*(0.42 + 0.0071*tree1->GetLeaf("ScalarEt")->GetValue(0));//*ratiotest->GetBinContent(ratiotest->FindBin(tree1->GetLeaf("MBScalarEt")->GetValue(0)));//1.4;
//cout<<pmcsmbz<<endl;
// pmcsmbz = pmcsmbz/(-0.00249 - (0.0003136*pmcsmbz) + (0.00121*pmcsmbz*pmcsmbz));
// if (pmcsmbz < 50) pmcsmbz = pmcsmbz/(-0.00005 - (0.00196*pmcsmbz) + (0.00132*pmcsmbz*pmcsmbz));
//cout<<pmcsmbz<<" new "<<endl;
//below way over corrects... obviously not compensating for "other" correctly...
// if (pmcsmbz > 0) pmcsmbz = pmcsmbz* ((0.5 * (1. + TMath::Erf((pmcsmbz-8.36)/(TMath::Sqrt(2)*4.9)))));
// if (pmcsmbz > 0) pmcsmbz = TMath::Sqrt(pmcsmbz);// * (1/(0.5 * (1. + TMath::Erf((pmcsmbz-8.36)/(sqrt(2)*4.9)))));
// else pmcsmbz = 0;
pmcsset = tree1->GetLeaf("ScalarEt")->GetValue(0) - tree1->GetLeaf("MBScalarEt")->GetValue(0) + pmcsmbz;
tpmcs->Fill(pmcsset, wgt);//minus mb plus mb
tpmcsd->Fill(pmcsset);
tpmcsother->Fill(pmcsset -pmcsmbz, wgt);
tpmcsmb->Fill(pmcsmbz, wgt);
pmcsmb->Fill(pmcsmbz, wgt);
tpmcssum->Fill(pmcsset, pmcsset);//minus mb plus mb
tpmcsothersum->Fill(pmcsset, pmcsset -pmcsmbz);
tpmcsmbsum->Fill(pmcsset, pmcsmbz);
test += pmcsset;
test2 += pmcsset -pmcsmbz;
test3 += pmcsmbz;
//cout<<pmcsset<<" "<<pmcsset -pmcsmbz<<" "<<pmcsmbz<<endl;
b->Fill(pmcsset, pmcsmbz);
t->Fill(pmcsset, pmcsmbz - pmcsset);
mbandset->Fill(pmcsset, pmcsmbz, wgt);
otherandset->Fill(pmcsset, pmcsset -pmcsmbz, wgt);
// tpmcs->Fill(pmcsmbz);
// tpmcs->Fill((pmcsset - pmcsmbz) + (pmcsmbz*(1.25 + 0.005*pmcsmbz- 0.0002*pmcsmbz*pmcsmbz)));
// tpmcs->Fill((pmcsset - pmcsmbz) + (pmcsmbz + (0.052 - 0.019*pmcsmbz+ 0.0018*pmcsmbz*pmcsmbz)));
// tpmcs->Fill((pmcsset - pmcsmbz) + (pmcsmbz * (0.052 - 0.019*pmcsmbz+ 0.0018*pmcsmbz*pmcsmbz)));
// tpmcs->Fill((pmcsset - pmcsmbz) + (1.6*pmcsmbz + (0.2 - 0.065*pmcsmbz+ 0.0023*pmcsmbz*pmcsmbz)));
//uncomment top before using
// tpmcs->Fill((pmcsset - pmcsmbz) + (pmcsmbz * (1.018 + 0.0319*pmcsmbz - 0.00228*pmcsmbz*pmcsmbz + 0.00008*pmcsmbz*pmcsmbz*pmcsmbz)));
}
}示例2: TrackTrackSignal
TH2D * TrackTrackSignal(double pttriglow , double pttrighigh , double ptasslow , double ptasshigh, int centmin, int centmax, double vzrange)
{
Long64_t nentries = c->GetEntries();
Long64_t nbytes = 0, nb = 0;
TH2D * hTrackTrackSignal = new TH2D(Form("signal_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#Delta#eta;#Delta#phi",100,-5,5,75,-pi,2*pi);
double ntottrig = 0;
int n_entries_in_cent_range = cent_index_start[centmax] - cent_index_start[centmin];
for (Long64_t jentry=cent_index_start[centmin]; jentry<cent_index_start[centmax];jentry++) {
if(jentry%1000==0) cout<<jentry-cent_index_start[centmin]<<"/"<<n_entries_in_cent_range<<endl;
c->GetEntry(jentry);
// if(c->evt->hiBin < centmin) continue;
// if(c->evt->hiBin >= centmax) break;
if(fabs(c->evt->vz)>vzrange) continue;
int ntrig = 0 , nass = 0;
vector<double> trigtrkEta;
vector<double> trigtrkPhi;
vector<double> trigtrkIndex;
vector<double> asstrkEta;
vector<double> asstrkPhi;
vector<double> asstrkIndex;
for(int i = 0 ; i < c->track->nTrk ; ++i)
{
if( c->track->trkPt[i]<pttriglow || c->track->trkPt[i]>pttrighigh || fabs(c->track->trkEta[i])>maxetatrg ) continue;
++ntrig;
trigtrkEta.push_back(c->track->trkEta[i]);
trigtrkPhi.push_back(c->track->trkPhi[i]);
trigtrkIndex.push_back(i);
}
for(int j = 0 ; j < c->track->nTrk ; ++j)
{
if( c->track->trkPt[j]<ptasslow || c->track->trkPt[j]>ptasshigh || fabs(c->track->trkEta[j])>maxetaass) continue;
++nass;
asstrkEta.push_back(c->track->trkEta[j]);
asstrkPhi.push_back(c->track->trkPhi[j]);
asstrkIndex.push_back(j);
}
for(int i = 0 ; i < ntrig ; ++i)
{
ntottrig += 1;
for(int j = 0 ; j < nass ; ++j)
{
if(trigtrkIndex[i]==asstrkIndex[j]) continue;
double deta = fabs(trigtrkEta[i]-asstrkEta[j]);
double dphi = fabs(trigtrkPhi[i]-asstrkPhi[j]);
if( dphi > pi ) dphi = 2*pi - dphi;
hTrackTrackSignal->Fill(deta,dphi);
hTrackTrackSignal->Fill(-deta,dphi);
hTrackTrackSignal->Fill(deta,-dphi);
hTrackTrackSignal->Fill(-deta,-dphi);
hTrackTrackSignal->Fill(deta,(2*pi)-dphi);
hTrackTrackSignal->Fill(-deta,(2*pi)-dphi);
}
}
// if(jentry>100) break;
}
hTrackTrackSignal->Scale(1/ntottrig);
return hTrackTrackSignal;
}示例3: Loop
void CaloHitAna::Loop(bool useProtonCut)
{
if (fChain == 0) return;
TFile *fcut = TFile::Open("qpid_cuts.root", "read");
if(!fcut)
{
cerr << "Could not open qpid_cuts.root for reading!\n";
return;
}
TCutG *cut_proton = dynamic_cast<TCutG*>(fcut->Get("qpid_proton"));
if(!cut_proton)
{
cerr << "Could not find TCutG qpid_proton!\n";
return;
}
TCutG *cut_pspx = dynamic_cast<TCutG*>(fcut->Get("pspx_zminus2"));
if(!cut_pspx)
{
cerr << "Could not find TCutG pspx_zminus1\n";
return;
}
fcut->Close();
delete fcut;
TH2D *hee = new TH2D("evse", "E1 vs E2", 200, 0, 400000, 200, 0, 400000);
TH2D *htt = new TH2D("thetavstheta", "Theta 1 vs Theta 2", 16, 27.2, 61.9, 16, 27.2, 61.9);
TH2D *hpp = new TH2D("phivsphi", "Phi 1 vs Phi 2", 180, 0, 180, 180, 180, 360);
TH1D *hsumtheta = new TH1D("hSumTheta", "Sum Theta", 180, 0, 180);
TH1D *hdiffphi = new TH1D("hDiffPhi", "Dela Phi", 180, 0, 180);
TH1D *hegamma = new TH1D("hEgamma", "Gamma Energy", 1000, 0, 10000);
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
int i;
int idxEmax1, idxEmax2;
Double_t sumTheta, eGamma;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
if(CaloHit_ < 2)
continue;
// Only take Z-1 events
if(!cut_pspx->IsInside(Pspx04e, Pspx05e))
continue;
idxEmax1 = -1;
idxEmax2 = -1;
for(i = 0; i < CaloHit_; i++)
{
// Only consider completely stopped protons
if(cut_proton->IsInside(CaloHit_fNf[i], CaloHit_fNs[i]) && CaloHit_fPhi[i] < TMath::Pi() && (idxEmax1 == -1 || CaloHit_fEnergy[i] > CaloHit_fEnergy[idxEmax1]))
{
// Petal 1
idxEmax1 = i;
}
else if(cut_proton->IsInside(CaloHit_fNf[i], CaloHit_fNs[i]) && CaloHit_fPhi[i] >= TMath::Pi() && (idxEmax2 == -1 || CaloHit_fEnergy[i] > CaloHit_fEnergy[idxEmax2]))
{
// Petal 2
idxEmax2 = i;
}
// if(idxEmax1 == -1 || CaloHit_fEnergy[i] > CaloHit_fEnergy[idxEmax1])
// {
// idxEmax2 = idxEmax1;
// idxEmax1 = i;
// }
// else if(idxEmax2 == -1 || CaloHit_fEnergy[i] > CaloHit_fEnergy[idxEmax2])
// idxEmax2 = i;
}
if(idxEmax1 != -1 && idxEmax2 != -1)
{
sumTheta = (CaloHit_fTheta[idxEmax1] + CaloHit_fTheta[idxEmax2]) * 180.0 / TMath::Pi();
htt->Fill(CaloHit_fTheta[idxEmax1] * 180.0 / TMath::Pi(), CaloHit_fTheta[idxEmax2] * 180.0 / TMath::Pi());
hpp->Fill(CaloHit_fPhi[idxEmax1] * 180.0 / TMath::Pi(), CaloHit_fPhi[idxEmax2] * 180.0 / TMath::Pi());
hsumtheta->Fill(sumTheta);
hdiffphi->Fill((CaloHit_fPhi[idxEmax2] - CaloHit_fPhi[idxEmax1]) * 180.0 / TMath::Pi());
if(sumTheta >= 70.0 && sumTheta <= 90.0)
{
hee->Fill(CaloHit_fEnergy[idxEmax1], CaloHit_fEnergy[idxEmax2]);
// Try to find coincident gammas
if(CaloHit_ > 2)
{
for(i = 0; i < CaloHit_; i++)
{
if(i != idxEmax1 && i != idxEmax2 && (eGamma = Doppler(CaloHit_fEnergy[i], CaloHit_fTheta[i])) < 10000)
// ToDo: Doppler correction
hegamma->Fill(eGamma);
}
}
}
//.........这里部分代码省略.........
示例4: Parity_pi
//.........这里部分代码省略.........
}
if(Day3<=138024) {
for(int jj=0; jj<Nrun_MB1; jj++) if(Day3 == bad_Ref_day3_MB1[jj]) {
Bad = 1;
break;
}
}
else if(Day3<=145020) {
for(int jj=0; jj<Nrun_MB2; jj++) if(Day3 == bad_Ref_day3_MB2[jj]) {
Bad = 1;
break;
}
}
else if(Day3<=154021) {}
else if(Day3<=165031) {
for(int jj=0; jj<Nrun_MB5; jj++) if(Day3 == bad_Ref_day3_MB5[jj]) {
Bad = 1;
break;
}
}
else {
for(int jj=0; jj<Nrun_MB6; jj++) if(Day3 == bad_Ref_day3_MB6[jj]) {
Bad = 1;
break;
}
}
if(Bad) continue; //bad run
if(RefMult) {
float gM = 0.9995 + 21.89/(4.191*RefMult-18.17) - 2.723e-5*(4.191*RefMult-18.17);
Eweight = gM + 0.0009326*(gM-1)*PVtxz;
}
hBz->Fill(Bz);
hTrigger->Fill(Trigger);
if((Trigger%5)<2) continue;
Centrality = 0;
for(int j=0; j<9; j++) if(RefMult>cenDef[j]) Centrality = j+1;
refmultCorrUtil.init(Run);
if ( refmultCorrUtil.isBadRun(Run) ) continue;
refmultCorrUtil.initEvent(RefMult, PVtxz, ZDCcoin) ;
Int_t cent9 = 1 + refmultCorrUtil.getCentralityBin9() ;
Eweight = refmultCorrUtil.getWeight();
if((i+1)%1000==0) cout<<Centrality<<" "<<Eweight<<endl;
Centrality = cent9;
hVertexZ->Fill(PVtxz);
hEventTally->Fill("Total Event",1);
if(TOFMult < -100+3.0*RefMult) continue; //remove pile-up
if(TOFMult > 180 +5.2*RefMult) continue;
if(TMath::Abs(PVtxz) > Vz_cut) continue; //Z-vertex cut; track quality cut done in PkNtupleMaker
if((PVtxz-VPDvz)>3 || (PVtxz-VPDvz)<-3) continue;
Ref_Day3->Fill(Day3,RefMult);
TOF_Day3->Fill(Day3,TOFMult);
NPT_Day3->Fill(Day3,NPTracks);
Hist_RefMult->Fill(RefMult);
Hist_TOFMult->Fill(TOFMult);
if(RefMult>10) {
p_RefMult->Fill(RefMult,RefMult,Eweight);
p_TOFMult->Fill(TOFMult,TOFMult,Eweight);
}
示例5: if
//.........这里部分代码省略.........
// cout<<mccommand<<endl;
if(mccommand==0)
{
if(!c->selectEvent()) continue;
if(!c->hlt->HLT_HIJet80_v1) continue;
}
if( fabs(c->evt->vz) > vzrange ) continue;
// if( c->akPu3PF->nref < 2 ) continue;
// double dijetdphi = fabs(c->akPu3PF->jtphi[0] - c->akPu3PF->jtphi[1]);
// if( dijetdphi > pi ) dijetdphi = 2*pi - dijetdphi;
// if( dijetdphi > dijet_dphi_CUT ) continue;
int jetindex = -1;
while(true)
{
// cout<<"here"<<endl;
jetindex++;
if( jetindex == c->akPu3PF->nref ) break;
if( c->akPu3PF->jtpt[jetindex] > leadingjetpthigh ) continue;
if( c->akPu3PF->jtpt[jetindex] < leadingjetptlow ) break;
if( fabs(c->akPu3PF->jteta[jetindex]) > jetamax || fabs(c->akPu3PF->jteta[jetindex]) < jetamin ) continue;
if ((c->akPu3PF->trackMax[jetindex]/c->akPu3PF->jtpt[jetindex])<0.01) continue;
// if(jentry%1000==0) cout<<"here"<<endl;
double jeteta = c->akPu3PF->jteta[jetindex];
double jetphi = c->akPu3PF->jtphi[jetindex];
int ntrig = 1 , nass = 0;
vector<double> asstrkEta(c->track->nTrk);
vector<double> asstrkPhi(c->track->nTrk);
vector<double> asstrkPt(c->track->nTrk);
vector<double> asstrkIndex(c->track->nTrk);
if(mccommand<2)
{
for(int j = 0 ; j < c->track->nTrk ; ++j)
{
if( c->track->trkPt[j]<ptasslow || c->track->trkPt[j]>ptasshigh || fabs(c->track->trkEta[j])>maxetaass) continue;
if( c->track->trkAlgo[j]<4 || c->track->highPurity[j])
{
++nass;
asstrkEta.push_back(c->track->trkEta[j]);
asstrkPhi.push_back(c->track->trkPhi[j]);
asstrkPt.push_back(c->track->trkPt[j]);
asstrkIndex.push_back(j);
}
}
}
else
{
for(int j = 0 ; j < c->track->nParticle ; ++j)
{
if( c->track->pPt[j]<ptasslow || c->track->pPt[j]>ptasshigh || fabs(c->track->pEta[j])>maxetaass) continue;
++nass;
asstrkEta.push_back(c->track->pEta[j]);
asstrkPhi.push_back(c->track->pPhi[j]);
asstrkPt.push_back(c->track->pPt[j]);
asstrkIndex.push_back(j);
}
}
for(int i = 0 ; i < ntrig ; ++i)
{
ntottrig += 1;
for(int j = 0 ; j < nass ; ++j)
{
double deta = fabs(jeteta-asstrkEta[j]);
double dphi = fabs(jetphi-asstrkPhi[j]);
// double ptw = asstrkPt[j];
if( dphi > pi ) dphi = 2*pi - dphi;
double effweight = 1;
if(doptweight!=0) effweight = c->getTrackCorrection(asstrkIndex[j]);
// if(jentry%1000==0) cout<<effweight<<endl;
// {
hJetTrackSignal->Fill(deta,dphi,effweight);
hJetTrackSignal->Fill(-deta,dphi,effweight);
hJetTrackSignal->Fill(deta,-dphi,effweight);
hJetTrackSignal->Fill(-deta,-dphi,effweight);
hJetTrackSignal->Fill(deta,(2*pi)-dphi,effweight);
hJetTrackSignal->Fill(-deta,(2*pi)-dphi,effweight);
// }
// else
// {
// hJetTrackSignal->Fill(deta,dphi,effweight);
// hJetTrackSignal->Fill(-deta,dphi,effweight);
// hJetTrackSignal->Fill(deta,-dphi,effweight);
// hJetTrackSignal->Fill(-deta,-dphi,effweight);
// hJetTrackSignal->Fill(deta,(2*pi)-dphi,effweight);
// hJetTrackSignal->Fill(-deta,(2*pi)-dphi,effweight);
// }
}
}
}
// if(jentry>100) break;
}
if( ntottrig != 0 ) hJetTrackSignal->Scale(1/ntottrig);
c->ResetBooleans();
return hJetTrackSignal;
}示例6: TCanvas
TCanvas *plot_occupancy(Int_t adc_cut=40, Int_t adc_neighbor_cut=10000, Int_t multiplicity_cut=12, Int_t tdc_min=750, Int_t tdc_width=300){
TString cut, draw, draw1, title;
title.Form("run_%d_Occupancy",run);
Int_t nbin=196;
Int_t min=1, max=197;
Int_t nbinm=11;
Int_t minm=-1, maxm=10;
TH1D *hoccupancy = new TH1D("hoccupancy","hoccupancy",nbin,min,max);
TH1D *hmultiplicity = new TH1D("hmultiplicity","hmultiplicity",nbinm,minm,maxm);
TH2D *heatmap = new TH2D("heatmap","heatmap",nbin,min,max,nbinm,minm,maxm);
TCanvas *cOCCUPANCY= new TCanvas("cOCCUPANCY",title,xcanvas,ycanvas);
Int_t nentries=n_events_to_analyze;
TString tmpentry;
MyStyle->SetStatX(0.9);
MyStyle->SetStatY(0.9);
MyStyle->SetStatW(0.2);
setPaddleIndices();
for (Int_t id=1;id<=nentries;id++){
T->GetEntry(id);
Int_t nmultiplicity=0;
Int_t good_paddle[500];
for(Int_t icount=0;icount<500;icount++){good_paddle[icount]=-1;}
for(Int_t pmt=0; pmt<NUMPMT;pmt++){
Int_t ipaddle = (pmt+1)*NUMPADDLE+1;
for(Int_t pixel=0; pixel<NUMPIXEL;pixel++){
Int_t index = pmt*NUMPIXEL+pixel;
if (pixel!=pixel1[pmt]-1&&pixel!=pixel2[pmt]-1){
ipaddle--;
if (tdcl[index]>tdc_min&&tdcl[index]<tdc_min+tdc_width){
nmultiplicity++;
if (ipaddle == 1){
if (adc_c[index]>adc_cut && adc_c[paddleindex[ipaddle+1]] < adc_neighbor_cut ){
good_paddle[nmultiplicity-1]=ipaddle;
}
}else if (ipaddle == NUMPMT*NUMPADDLE){
if (adc_c[index]>adc_cut && adc_c[paddleindex[ipaddle-1]] < adc_neighbor_cut ){
good_paddle[nmultiplicity-1]=ipaddle;
}
}else{
if (adc_c[index]>adc_cut && adc_c[paddleindex[ipaddle-1]] < adc_neighbor_cut && adc_c[paddleindex[ipaddle+1]] < adc_neighbor_cut ){
good_paddle[nmultiplicity-1]=ipaddle;
}
}
}
}
}
}
if(nmultiplicity>0&&nmultiplicity<=multiplicity_cut){
for(Int_t icount=0;icount<nmultiplicity;icount++){
hoccupancy->Fill(good_paddle[icount]);
heatmap->Fill(good_paddle[icount],nmultiplicity);
}
}
hmultiplicity->Fill(nmultiplicity);
}
cOCCUPANCY->Clear();
cOCCUPANCY->Divide(1,2) ;
title.Form("run_%d_OCCUPANCY_tdc_min_%d_max_%d.png",
run,tdc_min,tdc_min+tdc_width);
cOCCUPANCY->Print(title);
cOCCUPANCY->cd(1);
//gPad->SetLogy();
hoccupancy->Draw();
hoccupancy->GetXaxis()->SetNdivisions(NUMPMT,NUMPMT,0,0);
hoccupancy->SetLineColor(kBlue);
gPad->SetGridx();
TVirtualPad *c1_1 = cOCCUPANCY->cd(2);
c1_1->Divide(2,1);
c1_1->cd(1);
hmultiplicity->Draw();
hmultiplicity->SetLineColor(kBlue);
c1_1->cd(2);
heatmap->Draw("COLZ");
heatmap->GetXaxis()->SetNdivisions(NUMPMT,NUMPMT,0,0);
gPad->SetGridx();
return cOCCUPANCY;
}示例7: Loop
// LOOP function
// Note that looping through millions of events will take several minutes regardless of what you ask it to do with them.
void Loop(){
// These booleans are used to decide whether to add an EVENT (rather than a hit) to a histogram
bool left_showermax = false;
bool hit_showermax = false;
bool hit_lead = false;
bool left_lead = false;
// Here we define doubles to store data in
double primary_e_1 = 0;
double primary_r_1 = 0;
double exit_count_showermax = 0;
double exit_count_lead = 0;
double incident_count = 0;
double sum_showermax_energy = 0;
double sum_energy_post_lead = 0;
double sum_incident_energy = 0;
// basic check to make sure we actually pulled in something
if (fChain == 0) return;
// Counter used to print progress to terminal while the function is running
int events_crunched = 0;
// create a stopwatch, which starts timing
clock_t start_time = clock();
// This for loop iterates through all events, events can register hits but are not themselves hits
for (Long64_t event = 0; event < fChain->GetEntries(); fChain->GetEntry(event++)) {
// This is the timer code, which is extremely useful because it tells you relative information on how long it takes to process each event, and how changing geometry affects this.
// Furthermore it prevents the terminal from freezing for an hour, and lets you know quickly if the script is actually processing or if it's hung up on something.
events_crunched++;
int remainder = events_crunched % 100000;
if( remainder == 0){
clock_t runtime = clock();
clock_t net_time = (runtime - start_time);
double real_time = ((double)net_time)/CLOCKS_PER_SEC;
double avg_per_sec = events_crunched / real_time;
std::cout << " " << events_crunched << " events processed at average rate of " << avg_per_sec << " events per second." << std::endl;
};
// reset booleans
left_showermax = false;
hit_showermax = false;
primary_e_1 = 0;
primary_r_1 = 0;
exit_count_showermax = 0;
hit_lead = false;
left_lead = false;
incident_count = 0;
sum_showermax_energy = 0;
sum_energy_post_lead = 0;
sum_incident_energy = 0;
exit_count_lead = 0;
// Use 'if' statements to select what cuts you want on events, perform any analysis
// outside the loop in the function below.
//
// Looping should be linear time, passing through each event and hit only once.
// Script runtime will increase with each comparison
if( hit->size() > 0){
for (size_t i = 0; i < hit->size(); i++) {
// Ideally nest as many statements as possible to
// reduce absolute number of comparisons
//
// Note that Fill can take a second argument (or third for TH2) that is a weight. It seems like ev is defined as a remoll event in remolltypes.hh, so it's actually a struct with various attributes.
// The first cut should be the particle id in this case electrons, nest all other cuts for that particle within it to save time.
//
// *NOTE* positrons should be included in future simulations, simply change the line below to if(hit->at(i).pid == 11 || hit->at(i).pid == -11){
//
if(hit->at(i).pid == 11){
// Unless you intend to look at backsplashing particles, put a cut to only count particles moving forward in the lab frame i.e. pz > 0
// Front of showermax
if(hit->at(i).det == 30 && hit->at(i).pz > 0){
h_hit_e_showermax_front->Fill(hit->at(i).p, rate);
// A particle with track id 1 or 2 seems to refer to one of the moller scattered electrons, there are no particles with trid 0, 3, 4, or 5
// and trid 1 and 2 seem to correspond to mother track id 0 which means they came from the particle gun. This also places a 1 MeV minimum energy constraint
if((hit->at(i).trid == 1 || hit->at(i).trid == 2) && hit->at(i).p > 1){
// The precise z coordinates were implemented due to a bug, they can probably be removed without consequence since the bug was due to a very strange batch of data
// rather than anything in the actual code.
if(hit->at(i).z > 29073.98 && hit->at(i).z < 29074.02){
hit_showermax = true;
primary_e_1 = hit->at(i).p;
primary_r_1 = hit->at(i).r;
}
//.........这里部分代码省略.........
示例8: if
TH2D * HFTrackSignal(double pttriglow , double pttrighigh , double ptasslow , double ptasshigh, int centmin, int centmax, int nmin, int nmax, double vzrange)
{
Long64_t nentries = c->GetEntries();
Long64_t nbytes = 0, nb = 0;
TH2D * hTrackTrackSignal = new TH2D(Form("signal_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#Delta#eta;#Delta#phi",33,0,8,48,-pi,2*pi);
hmult = new TH1D(Form("hmult_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";N",300,0,300);
hpttrg = new TH1D(Form("hpttrg_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";p_{T}",100,0,10);
hetatrg = new TH1D(Form("hetatrg_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#eta",100,-5,5);
hphitrg = new TH1D(Form("hphitrg_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#phi",100,-pi,pi);
hmulttrg = new TH1D(Form("hmulttrg_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#N in p_{T} range",300,0,300);
hptass = new TH1D(Form("hptass_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";p_{T}",100,0,10);
hetaass = new TH1D(Form("hetaass_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#eta",100,-5,5);
hphiass = new TH1D(Form("hphiass_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow,(int)pttrighigh,(int)ptasslow,(int)ptasshigh,centmin,centmax),";#phi",100,-pi,pi);
cout<<"nmin: "<<nmin<<" nmax: "<<nmax<<endl;
hdz = new TH1D("hdz_highpurity",Form(";dz %d #leq N < %d",nmin,nmax),1000,-100,100);
hdz_dzerr = new TH1D("hdz_dzerr_highpurity",Form(";dz/dzerr %d #leq N < %d",nmin,nmax),1000,-50,50);
hdxy = new TH1D("hdxy_highpurity",Form(";dxy %d #leq N < %d",nmin,nmax),1000,-50,50);
hdxy_dxyerr = new TH1D("hdxy_dxyerr_highpurity",Form(";dxy/dxyerr %d #leq N < %d",nmin,nmax),1000,-50,50);
hpterr_pt = new TH1D("hpterr_pt_highpurity",Form(";pterror/pt %d #leq N < %d",nmin,nmax),1000,0,1);
hnhits = new TH1D("hnhits_highpurity",Form(";nhits %d #leq N < %d",nmin,nmax),80,0,80);
hchi2ndof = new TH1D("hchi2ndof_highpurity",Form(";chi2/ndof %d #leq N < %d",nmin,nmax),1000,0,10);
hetaphi = new TH2D("hetaphi_highpurity",Form(";#eta %d #leq N < %d;#phi",nmin,nmax),50,-2.5,2.5,62,-pi,pi);
int nmult = 0;
int n_entries_in_cent_range = cent_index_start[centmax] - cent_index_start[centmin];
// for (Long64_t jentry=cent_index_start[centmin]; jentry<cent_index_start[centmax];jentry++) {
for (Long64_t jentry=0; jentry<nentries;jentry++) {
if(jentry%1000==0)
{
cout<<jentry<<"/"<<nentries<<endl;
}
// if(jentry==9999) break;
c->GetEntry(jentry);
// if(c->evt.hiBin < centmin) continue;
// if(c->evt.hiBin >= centmax) break;
// if(!(c->skim.phfPosFilter1&&c->skim.phfNegFilter1&&c->skim.phltPixelClusterShapeFilter&&c->skim.pprimaryvertexFilter)) continue;
if(skipevent(vzrange,211256)) continue;
// if(c->evt.run!=202792) continue;
// if(c->evt.lumi<317 || c->evt.lumi>1014) continue;
int thismult = 0;
for(int i = 0 ; i < c->track.nTrk ; ++i)
{
if(c->track.trkPt[i]>0.4&&fabs(c->track.trkEta[i])<2.4&&c->track.highPurity[i]&&fabs(c->track.trkDz1[i]/c->track.trkDzError1[i])<3&&fabs(c->track.trkDxy1[i]/c->track.trkDxyError1[i])<3&&c->track.trkPtError[i]/c->track.trkPt[i]<0.1) thismult++;
}
hmult->Fill(thismult);
// cout<<thismult<<endl;
if(thismult<nmin || thismult>=nmax)
{
cout<<"WARNING: different mult in input file. nmin "<<nmin<<", nmax "<<nmax<<", thismult "<<thismult<<", nTrk "<<c->track.nTrk<<endl;
continue;
}
nmult++;
int ntrig = 0 , nass = 0;
vector<double> trigtrkEta;
vector<double> trigtrkPhi;
vector<double> trigtrkIndex;
vector<double> asstrkEta;
vector<double> asstrkPhi;
vector<double> asstrkIndex;
for(int i = 0 ; i < c->track.nTrk ; ++i)
{
hetaphi->Fill(c->track.trkEta[i],c->track.trkPhi[i]);
hdz->Fill(c->track.trkDz1[i]);
hdxy->Fill(c->track.trkDxy1[i]);
hdz_dzerr->Fill(c->track.trkDz1[i]/c->track.trkDzError1[i]);
hdxy_dxyerr->Fill(c->track.trkDxy1[i]/c->track.trkDxyError1[i]);
hpterr_pt->Fill(c->track.trkPtError[i]/c->track.trkPt[i]);
hnhits->Fill(c->track.trkNHit[i]);
hchi2ndof->Fill(c->track.trkChi2[i]/c->track.trkNdof[i]);
if( c->track.trkPt[i]<pttriglow || c->track.trkPt[i]>pttrighigh || fabs(c->track.trkEta[i])>maxetatrg ) continue;
bool considertrack = false;
if(mytrackquality==0&&fabs(c->track.trkEta[i])<2.4&&c->track.highPurity[i]&&fabs(c->track.trkDz1[i]/c->track.trkDzError1[i])<3&&fabs(c->track.trkDxy1[i]/c->track.trkDxyError1[i])<3&&c->track.trkPtError[i]/c->track.trkPt[i]<0.1)
{
considertrack=true;
}
else if(mytrackquality==1)
{
considertrack=true;
}
else if(mytrackquality==2&&!(fabs(c->track.trkEta[i])<2.4&&c->track.highPurity[i]&&fabs(c->track.trkDz1[i]/c->track.trkDzError1[i])<3&&fabs(c->track.trkDxy1[i]/c->track.trkDxyError1[i])<3&&c->track.trkPtError[i]/c->track.trkPt[i]<0.1))
{
considertrack=true;
}
if(considertrack)
{
++ntrig;
trigtrkEta.push_back(c->track.trkEta[i]);
trigtrkPhi.push_back(c->track.trkPhi[i]);
trigtrkIndex.push_back(i);
hpttrg->Fill(c->track.trkPt[i]);
//.........这里部分代码省略.........
示例9: draw5
void draw5(){
gStyle->SetOptStat("mre");
// TFile* f = TFile::Open("testing/corInterp/latestNums/generateAngularResolutionTreePlots_352_2016-02-21_15-38-05.root"); //generateAngularResolutionTreePlots_352_2016-02-19_18-01-52.root");
TChain* c = new TChain("angResTree");
UInt_t l3TrigPatternH;
c->SetBranchAddress("l3TrigPatternH", &l3TrigPatternH);
Double_t phiExpected;
c->SetBranchAddress("phiExpected", &phiExpected);
// c->Add("newLindaNumbers_4steps_WAISHPOL_NEW11_cosminV3_nfixedBug_2016_02_05_time_15_42_15/*.root");
c->Add("newLindaNumbers_4steps_WAISHPOL_NEW11_cosminV3_nfixedBug_2016_02_05_time_15_42_15/generateAngularResolutionTreePlots_352*.root");
const Long64_t numEntries = c->GetEntries();
ProgressBar p(numEntries);
TH2D* h = new TH2D("h", "#phi_{expected} vs. L3 triggered #phi-sectors; #phi_{expected} (Degrees); #phi-sector", 360, 0, 360, NUM_PHI, 0, NUM_PHI);
TH1D* h2 = new TH1D("h2", "Angle between triggered #phi-sector center and #phi_{expected}; #delta#phi (Degrees); Events / bin", 360, -180, 180);
TH1D* h3 = new TH1D("h3", "#phi-sectors triggered; #phi-sector; Number of triggers", NUM_PHI, 0, NUM_PHI);
TH1D* h4 = new TH1D("h4", "#delta#phi-sectors; #phi-sector; Number of triggers", 5, -2, 3);
for(Long64_t entry=0; entry < numEntries; entry++){
c->GetEntry(entry);
for(int phi=0; phi<NUM_PHI; phi++){
Int_t l3 = RootTools::getBit(phi, l3TrigPatternH);
if(l3){
h->Fill(phiExpected, phi);
Double_t deltaPhiThingy = RootTools::getDeltaAngleDeg(phiExpected + 45, 22.5*phi);
h2->Fill(deltaPhiThingy);
h3->Fill(phi);
if(deltaPhiThingy <= -11.25){
h4->Fill(-1);
}
else if(deltaPhiThingy > 11.25){
h4->Fill(1);
}
else{
h4->Fill(0);
}
}
}
p++;
}
new TCanvas();
h->Draw("colz");
new TCanvas();
h2->Draw();
new TCanvas();
h4->Draw();
for(int binx=1; binx<=h4->GetNbinsX(); binx++){
h4->GetXaxis()->SetBinLabel(binx, TString::Format("%d", -3 + binx));
}
h4->GetXaxis()->SetLabelOffset(.01);
}示例10: GradeCorrelation
//.........这里部分代码省略.........
// Find next regular term only
bool foundRegTerm = false;
for (int jTerm = iTerm + 1; jTerm < nTerms && !foundRegTerm; ++jTerm) {
const Student::Enrollment jEnrollment = student->Enrollments()[jTerm];
if (!MyFunctions::regularSemester(jEnrollment.term)) continue;
foundRegTerm = true;
for (Student::Grade jGrade : jEnrollment.grades) {
if (jGrade.course == iGrade.course) continue;
if (!MyFunctions::ValidGrade(jGrade.grade)) continue;
double prediction_j = student->CourseGradePrediction(jGrade, Student::DISTRIBUTION);
double delta_j = jGrade.quality - prediction_j;
corrMap[std::make_pair(iGrade.course, jGrade.course)].Add(delta_i, delta_j);
++nPairAll;
}
}
}
}
}
myBenchmark->Stop("Main Loop");
std::cout << "nPairAll = " << nPairAll << std::endl;
std::cout << "Unique Pairs = " << corrMap.size() << std::endl;
TH1D* rHist = new TH1D("rHist", "Correlation Coefficient, #rho", 120, -1.2, 1.2);
TH1D* pHist = new TH1D("pHist", "Probablity Distribution", 100, 0., 1.);
TH1D* nHist = new TH1D("nHist", "Number of entries", 100, 0., 2000.);
TH2D* pVrHist = new TH2D("pVrHist", "Prob vs. #rho", 100, -1., 1., 100, 0., 1.);
myBenchmark->Start("Prune");
for (auto iter = corrMap.begin(); iter != corrMap.end();) {
if (iter->second.n() < nCut) {
corrMap.erase(iter++);
continue;
}
double p = iter->second.p();
double r = iter->second.r();
// Test for nan?
if (p != p) {
std::cout << "Found p = nan: n = " << iter->second.n() << std::endl;
corrMap.erase(iter++);
continue;
}
if (p < 0.) {
corrMap.erase(iter++);
continue;
}
rHist->Fill(r);
pHist->Fill(p);
nHist->Fill(iter->second.n());
pVrHist->Fill(r, p);
if (p < 1. - prob && p > prob) {
corrMap.erase(iter++);
}
else {
std::cout << "r = " << iter->second.r() << ", p = " << p << std::endl;
++iter;
}
}
myBenchmark->Stop("Prune");
std::cout << "Post Cut = " << corrMap.size() << std::endl;
myBenchmark->Start("Sort");
std::vector<std::pair<std::pair<TString, TString>, CorrelationCalculator>> corrVec(corrMap.begin(), corrMap.end());
std::sort(corrVec.begin(), corrVec.end(), &sortFunc);
myBenchmark->Stop("Sort");
int printTop = 50;
int printed = 0;
for (auto const& entry : corrVec) {
std::cout << entry.first.first << " : " << entry.first.second << "\t, n = " << entry.second.n() << "\t, r = " << entry.second.r()
<< "\t, p = " << entry.second.p() << std::endl;
++printed;
if (printed >= printTop) break;
}
TCanvas* c1 = new TCanvas("c1", "Grade Correlation", 1600, 1200);
c1->Divide(2,2);
c1->cd(1);
TF1* myGaus = new TF1("myGaus", "gaus", -1., 1.);
myGaus->SetParameters(600., 0., 0.2);
myGaus->FixParameter(1, 0.);
rHist->Fit(myGaus, "0B", "", -1., 0.);
rHist->DrawCopy();
myGaus->DrawCopy("SAME");
c1->cd(2);
pHist->DrawCopy();
c1->cd(3);
nHist->DrawCopy();
c1->cd(4);
pVrHist->DrawCopy();
float rt, cp;
myBenchmark->Summary(rt, cp);
}示例11: L1q1trigger
void L1q1trigger(){
const TString l1_input = "~/scratch1/DiHadronCorrelations/L1UpgradeAnalyzer.root";
TFile *lFile = TFile::Open(l1_input);
Int_t l1Up_evt, l1Up_run, l1Up_et, l1Up_q2;
TTree *l1UpTree = (TTree*)lFile->Get("L1UpgradeAnalyzer/L1UpgradeTree");
l1UpTree->SetBranchStatus("*",0);
l1UpTree->SetBranchStatus("event",1);
l1UpTree->SetBranchStatus("run",1);
l1UpTree->SetBranchStatus("centrality_hwPt",1);
l1UpTree->SetBranchStatus("v2_hwPt",1);
l1UpTree->SetBranchAddress("event",&l1Up_evt);
l1UpTree->SetBranchAddress("run",&l1Up_run);
l1UpTree->SetBranchAddress("centrality_hwPt",&l1Up_et);
l1UpTree->SetBranchAddress("v2_hwPt",&l1Up_q2);
const TString forest_input = "~/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/ebyeflow_promptreco_pt009.root";
TFile *fFile = TFile::Open(forest_input);
TTree *fTree = (TTree*)fFile->Get("ebyeflow_ana_HI/q2ntuple");
Float_t f_evt, f_run, f_lumi, hiBin;
float hiHF,q2hf;
fTree->SetBranchStatus("*",0);
fTree->SetBranchStatus("evtnum",1);
fTree->SetBranchStatus("run",1);
fTree->SetBranchStatus("lumi",1);
fTree->SetBranchStatus("cent",1);
fTree->SetBranchStatus("q2hf",1);
fTree->SetBranchStatus("hf",1);
fTree->SetBranchAddress("evtnum",&f_evt);
fTree->SetBranchAddress("run",&f_run);
fTree->SetBranchAddress("lumi",&f_lumi);
fTree->SetBranchAddress("cent",&hiBin);
fTree->SetBranchAddress("hf",&hiHF);
fTree->SetBranchAddress("q2hf",&q2hf);
map<long, int> kmap;
map<long, int> kmapcal;
cout << "Filling the map..." << endl;
int l1up_entries = l1UpTree->GetEntries();
//int l_entries = lEvtTree->GetEntries();
for(long j = 0; j < l1up_entries; ++j){
if(j % 10000 == 0) printf("%ld / %d\n",j,l1up_entries);
l1UpTree->GetEntry(j);
long key = makeKey(l1Up_run, l1Up_evt);
pair<long,int> p(key,j);
kmap.insert(p);
kmapcal.insert(p);
}
cout << "map filled!!!" << endl;
//q2 histos
TH2D* q2CorrHist = new TH2D("q2CorrHist",";offline q_{2};online q_{2}",1000,0.,1.0,1000,0.0,1.0);
TH2D* HFCorrHist = new TH2D("HFCorrHist",";offline HF E_{T} sum (GeV); online HF E_{T} sum (GeV)",100,0.,8000.0,100,0.0,8000.0);
TH2D* q2HFCorrHist = new TH2D("q2HFCorrHist",";online HF*q^{2}_{2} (GeV); online HF E_{T} sum (GeV)",300,0.,30000.0,100,0.0,8000.0);
Float_t q2On = 0.0;
int entries = fTree->GetEntries();
for(long int j = 1; j < entries; ++j){
if(j % 10000 == 0) printf("%ld / %d \n",j,entries);
fTree->GetEntry(j);
long keycal = makeKey(f_run, f_evt);
map<long,int>::const_iterator gotcal = kmapcal.find(keycal);
if(gotcal == kmapcal.end()){
continue;
}
else {
l1UpTree->GetEntry(gotcal->second);
kmapcal.erase(keycal);
if(l1Up_et == 0) continue;
//q2 part
q2On = sqrt(l1Up_q2)/l1Up_et;
// if(l1Up_run<182060 && l1Up_run>182050 & l1Up_et<2839 && l1Up_et>1198)
// if(l1Up_et<2839 && l1Up_et>1198)
{
q2CorrHist->Fill(q2hf,q2On);
HFCorrHist->Fill(hiHF,l1Up_et);
q2HFCorrHist->Fill(l1Up_q2,l1Up_et);
}
q2On = 0.0;
}
}
TCanvas* c = new TCanvas("c","c",800,400);
c->Divide(2,1);
c->cd(1);
c->GetPad(1)->SetLogz();
q2CorrHist->Rebin2D(2,2);
//.........这里部分代码省略.........
示例12: Loop
// LOOP function
// Note that looping through millions of events will take several minutes regardless of what you ask it to do with them.
void Loop(){
double primary_e_1 = 0;
double primary_r_1 = 0;
double exit_count_showermax = 0;
double exit_count_lead = 0;
double incident_count = 0;
double sum_showermax_energy = 0;
double sum_energy_post_lead = 0;
double sum_incident_energy = 0;
if (fChain == 0) return;
// Counter used to print progress to terminal while the function is running
int events_crunched = 0;
// create a stopwatch, which starts timing
clock_t start_time = clock();
// This for loop iterates through all events, events can register hits but are not themselves hits
for (Long64_t event = 0; event < fChain->GetEntries(); fChain->GetEntry(event++)) {
events_crunched++;
int remainder = events_crunched % 100000;
if( remainder == 0){
clock_t runtime = clock();
clock_t net_time = (runtime - start_time);
double real_time = ((double)net_time)/CLOCKS_PER_SEC;
double avg_per_sec = events_crunched / real_time;
std::cout << " " << events_crunched << " events processed at average rate of " << avg_per_sec << " events per second." << std::endl;
};
// reset doubles
primary_e_1 = 0;
primary_r_1 = 0;
exit_count_showermax = 0;
incident_count = 0;
sum_showermax_energy = 0;
sum_energy_post_lead = 0;
sum_incident_energy = 0;
exit_count_lead = 0;
// Use 'if' statements to select what cuts you want on events, perform any analysis
// outside the loop in the function below.
//
// Looping should be linear time, passing through each event and hit only once.
// Script runtime will increase with each comparison
if( hit->size() > 0){
for (size_t i = 0; i < hit->size(); i++) {
// Ideally nest as many statements as possible to
// reduce absolute number of comparisons
//
// Note that Fill can take a second argument (or third for TH2) that is a weight. It seems like ev is defined as a remoll event in remolltypes.hh, so it's actually a struct with various attributes.
// // Could include pi+'s however it's probably not even worth the extra check every loop.
if(hit->at(i).pid == -211){
// Main Detector
if(hit->at(i).det == 28){
if(hit->at(i).pz > 0 && hit->at(i).r > 600){
h_hit_r_moller_det->Fill(hit->at(i).r, rate);
};
}
// Front of showermax
else if(hit->at(i).det == 30 && hit->at(i).pz > 0){
h_hit_e_showermax_front->Fill(hit->at(i).p, rate);
// A particle with track id 1 or 2 seems to refer to one of the moller scattered electrons, there are no particles with trid 0, 3, 4, or 5
// and trid 1 and 2 seem to correspond to mother track id 0 which means they came from the particle gun. This also places a 1 MeV minimum energy constraint
if((hit->at(i).trid == 1 || hit->at(i).trid == 2) && hit->at(i).p > 1){
if(hit->at(i).z > 29073.98 && hit->at(i).z < 29074.02){
primary_e_1 = hit->at(i).p;
primary_r_1 = hit->at(i).r;
// Note that an x-y plot was added in this script
h_x_y_showermax_incident->Fill(hit->at(i).x,hit->at(i).y,rate);
}
}
}
// Back of showermax
else if(hit->at(i).det == 31 && hit->at(i).pz > 0 && hit->at(i).p > 1){
// because the booleans have been removed we need to fill the plots at their respective detectors, this would be the typical method for an analysis script
exit_count_showermax++;
sum_showermax_energy += hit->at(i).p;
h_ev_e_showermax_front->Fill(primary_e_1,rate);
h_ev_r_showermax_front->Fill(primary_r_1,rate);
}
// pre-lead
else if(hit->at(i).det == 4050 && hit->at(i).pz > 0 && hit->at(i).p > 1){
// Again, booleans have been removed so plots are filled as they come up
incident_count++
sum_incident_energy += hit->at(i).p;
h_x_y_lead_incident->Fill(hit->at(i).x,hit->at(i).y,rate);
h_ev_e_showermax_thru->Fill(primary_e_1,rate);
h_ev_r_showermax_thru->Fill(primary_r_1,rate);
h_ev_showermax_exit_from_primary_hits->Fill(exit_count_showermax,rate);
//.........这里部分代码省略.........
示例13: plot_xtalE
void plot_xtalE(){
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);
gStyle->SetTitleSize(0.08,"t");
// a list of runs corresponding to xtal0 to xtal53
int runlist[54] = {3295,3296,3297,3369,3300,3301,3302,3303,3368,
3244,3245,3246,3247,3248,3249,3250,3251,3252,
3240,3239,3238,3237,3236,3235,3232,3234,3233,
3267,3266,3265,3263,3258,3257,3256,3254,3253,
3268,3272,3274,3275,3276,3277,3278,3281,3283,
3294,3293,3292,3289,3288,3287,3286,3285,3284};
// for(int i=0;i<54;i++){
// runlist[i] = 3305;
// }
ofstream outfile("energycalibration.fcl");
TCanvas *c1 = new TCanvas("c1","c1",1800,1200);
TH2D *xtalEnergy = new TH2D("xtalEnergy","xtalEnergy",9,0,9,6,0,6);
outfile << "// constants from runs" <<endl;;
for(int i=0;i<54;i++){
outfile << runlist[i] << " ";
if(i%9==8 && i!=53) {
outfile << endl;
outfile << "// ";}
}
outfile << "\nBEGIN_PROLOG" << endl;
outfile << "xtal_energy_response: {" << endl;
double mean[54];
double rms[54];
double norm[54];
double index[54];
c1->Divide(9,6);
for(int i=0;i<54;i++){
index[i]=i;
c1->cd(54-i);
c1->cd(54-i)->SetLogz();
//if(i==3 || i==8 || i==9) continue; // these are runs where the beams are at the wrong place
TFile *file = new TFile(Form("run%d.root",runlist[i]));
TH2D *beamEnergy2D = (TH2D*)file->Get("beamEnergy2D");
beamEnergy2D->Draw("colz text");
beamEnergy2D->SetMarkerSize(1);
beamEnergy2D->SetMaximum(3000);
beamEnergy2D->SetMinimum(1);
beamEnergy2D->SetTitle(Form("Run %d",runlist[i]));
double maxE = beamEnergy2D->GetMaximum();
int maxBin = beamEnergy2D->GetMaximumBin();
int x,y,z;
beamEnergy2D->GetBinXYZ(maxBin,x,y,z);
cout<<x<<" "<<y<<" "<<z<<" "<<(y-1)*9+9-x<<endl;
int xtalNum = (y-1)*9+9-x;
TH1D *beamEnergy = (TH1D*)file->Get(Form("beamEnergy%02d",xtalNum));
TH1D *beamTime = (TH1D*)file->Get(Form("beamTime%02d",xtalNum));
TH1D *beamTimeEnergy = (TH1D*)file->Get(Form("beamTimeEnergy%02d",xtalNum));
TH1D *syncEnergy = (TH1D*)file->Get(Form("syncEnergy%02d",xtalNum));
beamEnergy->Draw();
TF1 *fit = new TF1("fit","gaus(0)",1500,2800);
fit->SetParLimits(1,1700,2500);
fit->SetParLimits(2,50,150);
fit->SetParameters(100,2100,100);
beamEnergy->Fit("fit","REM");
norm[i]=fit->GetParameter(0);
mean[i]=fit->GetParameter(1);
rms[i]=fit->GetParameter(2);
TF1 *refit = new TF1("refit","gaus(0)",mean[i]-2*rms[i], mean[i]+2*rms[i]);
refit->SetParameters(norm[i],mean[i],rms[i]);
beamEnergy->Fit("refit","REM");
norm[i]=refit->GetParameter(0);
mean[i]=refit->GetParameter(1);
rms[i]=refit->GetParameter(2);
outfile << Form("xtal%d : %f",i,mean[i]) << endl;
xtalEnergy->Fill(8-i%9,i/9,mean[i]);
TText *text = new TText(0.15,0.75,Form("E%d=%.1f",xtalNum,mean[i]));
text->SetTextSize(0.09);
text->SetTextColor(2);
text->SetNDC();
beamEnergy->GetListOfFunctions()->Add(text);
}
//.........这里部分代码省略.........
示例14: Eloss
void Eloss()
{
TString dir = getenv("VMCWORKDIR");
TString fFileNamebase;
std::ifstream* ElossData;
fFileNamebase = dir+"/macro/Simulation/data/109Pd_Eloss.txt";
ElossData = new std::ifstream(fFileNamebase.Data());
if(ElossData->fail()){
std::cout<<cRED<<" = No file found! Please, check the path. Current :"<<fFileNamebase.Data()<<cNORMAL<<std::endl;
}
TCanvas *c1 = new TCanvas();
c1->Draw();
c1->Divide(1,2);
TH2D *Bragg = new TH2D("Bragg", "Bragg", 1000, 0,500, 1000, 0, 1);
TH2D *Interpolate = new TH2D("Interpolate","Interpolate", 1000, 0,200000, 100, 0, 1);
Double_t Eloss_elec, Eloss_nuc,SRIM_energy;
Double_t x = 0; //microns
Double_t x1 = 0;
Double_t x0 = 0;
Double_t Z = 46;
Double_t A = 109;
Double_t step = 1; //microns
Double_t energy = 190000; //keV
Double_t e1, e2, e3, dE1, dE2, dE3;
Double_t param0, param1, param2,param3,param4,param5,param6;
std::vector<std::tuple<Double_t,Double_t,Double_t>> SRIM_table;
std::cout<<cGREEN<<"Extracting data. Please wait...."<<cNORMAL<<std::endl;
while (!ElossData->eof()){ //Fill vector with datafile
*ElossData>>SRIM_energy>>Eloss_elec>>Eloss_nuc; //keV >> keV/micron >> keV/micron
SRIM_table.push_back(std::make_tuple(SRIM_energy, Eloss_elec, Eloss_nuc));
Interpolate->Fill(SRIM_energy, Eloss_elec);
}//EndWhile
/*for (Int_t n = 0; n < SRIM_table.size(); n++){
std::cout <<' '<< std::get<0> (SRIM_table.at(n))<<' '<<std::get<1> (SRIM_table.at(n))<<' '<<std::get<2> (SRIM_table.at(n))<<std::endl;
}//EndFor
*/
/*Int_t counter = SRIM_table.size()-3;
Double_t nuc_Loss, elec_Loss, total_Loss;
e1 = std::get<0> (SRIM_table.at(counter));
e2 = std::get<0> (SRIM_table.at(counter-1));
e3 = std::get<0> (SRIM_table.at(counter-2));
dE1 = std::get<1> (SRIM_table.at(counter));
dE2 = std::get<1> (SRIM_table.at(counter-1));
dE3 = std::get<1> (SRIM_table.at(counter-2));
Interpolate->Fill(e1, dE1);
Interpolate->Fill(e2, dE2);
Interpolate->Fill(e3, dE3);
Interpolate->Fit("enFit");
param0 = enFit->GetParameter(0);
param1 = enFit->GetParameter(1);
while (energy>20){
if(energy>std::get<0> (SRIM_table.at(counter))){
elec_Loss = param0*energy+param1;
x+=step;
//nuc_Loss = (std::get<2> (SRIM_table.at(counter+1)))*step;
//elec_Loss = (std::get<1> (SRIM_table.at(counter+1)))*step;
total_Loss = elec_Loss; //+nuc_Loss;
energy -= total_Loss;
Bragg->Fill(x/1000, total_Loss*10);
std::cout<<energy<<" "<<total_Loss<<" "<<x<<std::endl;
}
else{
counter-=1;
Interpolate->Reset();
e1 = std::get<0> (SRIM_table.at(counter));
e2 = std::get<0> (SRIM_table.at(counter-1));
e3 = std::get<0> (SRIM_table.at(counter-2));
dE1 = std::get<1> (SRIM_table.at(counter));
dE2 = std::get<1> (SRIM_table.at(counter-1));
dE3 = std::get<1> (SRIM_table.at(counter-2));
Interpolate->Fill(e1, dE1);
Interpolate->Fill(e2, dE2);
Interpolate->Fill(e3, dE3);
Interpolate->Fit("enFit");
param0 = enFit->GetParameter(0);
param1 = enFit->GetParameter(1);
}
}//EndWhile
*/
TF1 *enFit = new TF1("enFit", "[0]+[1]*x+[2]*pow(x,2)+[3]*pow(x,3)+[4]*pow(x,4)+[5]*pow(x,5)+[6]*pow(x,6)", 0 , 200000);
enFit->SetParameters(1,1,1,1,1,1);
Interpolate->Fit("enFit");
param0 = enFit->GetParameter(0);
param1 = enFit->GetParameter(1);
param2 = enFit->GetParameter(2);
param3 = enFit->GetParameter(3);
param4 = enFit->GetParameter(4);
param5 = enFit->GetParameter(5);
param6 = enFit->GetParameter(6);
//.........这里部分代码省略.........
示例15: main
//.........这里部分代码省略.........
m1max = 1300;
m2step = 2;
m2min = 2;
m2max = 26;
} else if (string(argv[1]) == "bprime") {
limitFile = "/cms/thomassen/2012/EWKino/Statistics/bprime/limits.txt";
outFile = "exclusions_bprime.root";
m1step = 50;
m1min = 500;
m1max = 700;
m2step = 0.1;
m2min = 0;
m2max = 1;
} else if (string(argv[1]) == "StopHiggs") {
limitFile = "limits/slide_combined_full_noreweight_chargino250.out";
outFile = "exclusions/exclusions_slide_combined_full_noreweight_chargino250.root";
m1step = 50;
m1min = 200;
m1max = 800;
m2step = 0.1;
m2min = 0;
m2max = 1;
} else {
cout << "Unknown keyword " << argv[1] << endl;
return -1;
}
if(m1max < m1min || m2max < m2min) return -1;
int m1bins = (int)((m1max - m1min) / m1step + 1);
int m2bins = (int)((m2max - m2min) / m2step + 1);
m1max = m1min + (m1bins - 1) * m1step;
m2max = m2min + (m2bins - 1) * m2step;
TH2D* hObs = new TH2D ("hObs", limitFile.c_str()
, m1bins, m1min - m1step/2., m1max + m1step/2.
, m2bins, m2min - m2step/2., m2max + m2step/2.);
hObs->GetXaxis()->SetTitle("m_{1}");
hObs->GetYaxis()->SetTitle("m_{2}");
TH2D* hExp = new TH2D(*hObs);
hExp->SetName("hExp");
TH2D* hExp1p = new TH2D(*hObs);
hExp1p->SetName("hExp1p");
TH2D* hExp1m = new TH2D(*hObs);
hExp1m->SetName("hExp1m");
TH2D* hExp2p = new TH2D(*hObs);
hExp2p->SetName("hExp2p");
TH2D* hExp2m = new TH2D(*hObs);
hExp2m->SetName("hExp2m");
cout << "[PT] Reading limits from " << limitFile << " ..." << endl;
ifstream in (limitFile.c_str());
while(in) {
string line;
getline(in, line);
if(!in) break;
double m1;
double m2;
float limit[6];
if(line[0] == '#') continue;
sscanf (line.c_str(),"%lf:%lf cls: %f %f ( %f : %f ) ( %f : %f )", &m1, &m2, limit, limit+1, limit+2, limit+3, limit+4, limit+5);
//cout << m1 << " " << m2 << " " << limit[0] << " " << limit[1] << " " << limit[2] << " " << limit[3] << " " << limit[4] << " " << limit[5] << endl;
if(m1 < m1min || m2 < m2min || m1 > m1max || m2 > m2max) continue;
if(!(m1 >= m1min && m1 <= m1max && (fmod(m1 - m1min, m1step) == 0 || fmod(m1 - m1min, m1step) - m1step < 1E-12)) || !(m2 >= m2min && m2 <= m2max && (fmod(m2 - m2min, m2step) == 0 || fmod(m2 - m2min, m2step) - m2step < 1E-12))) {
cout << "[PT] Skipping point (" << m1 << "," << m2 << ") because it is not on the grid" << endl;
continue;
}
if(hObs->GetBinContent(hObs->FindFixBin(m1, m2)) != 0 || hExp->GetBinContent(hExp->FindFixBin(m1, m2)) != 0) {
cout << "[PT] ERROR: The point (" << m1 << "," << m2 << ") was already filled." << endl;
return -1;
}
hObs->Fill(m1, m2, limit[0]);
hExp->Fill(m1, m2, limit[1]);
hExp1p->Fill(m1, m2, limit[2]);
hExp1m->Fill(m1, m2, limit[3]);
hExp2p->Fill(m1, m2, limit[4]);
hExp2m->Fill(m1, m2, limit[5]);
}
TFile f(outFile.c_str(), "RECREATE");
hObs->Write();
hExp->Write();
hExp1p->Write();
hExp1m->Write();
hExp2p->Write();
hExp2m->Write();
f.Close();
cout << "Output file produced: " << outFile << endl;
return 0;
}