本文整理汇总了C++中TH1I::Fill方法的典型用法代码示例。如果您正苦于以下问题:C++ TH1I::Fill方法的具体用法?C++ TH1I::Fill怎么用?C++ TH1I::Fill使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TH1I
的用法示例。
在下文中一共展示了TH1I::Fill方法的11个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: TCanvas
void A1::PlotChannel( unsigned aModule, unsigned aChannel ){
if (fChain == 0) return;
if( c1 ) delete c1;
c1 = new TCanvas( 3 );
stringstream stmp;
stmp << "Module " << aModule + 1 << ", Channel #" << aChannel << ";Time [0.1 ns]; Counts";
TH1I * hChannel = new TH1I( "hChannel", stmp.str().c_str(), 100000, 10000, 100000 );
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry);
nbytes += nb;
// if (Cut(ientry) < 0) continue;
vector<unsigned int> * HITS_PER_CHANNEL = HitsPerChannel[ aModule ];
vector<unsigned int> * CH = _Data[ aModule ][ aChannel ];//Module2_LE_CH32;
//std::cout << "HITS_PER_CHANNEL 32: " << HITS_PER_CHANNEL->at( 32 ) << std::endl;
if( HITS_PER_CHANNEL->at( aChannel ) > 0 ){
long trig = CH->at( 0 );
hChannel->Fill( trig );
}
}
hChannel->Draw();
}
示例2: strcpy
plotClasses(const char* canvas, const char* title,
int nClasses, const char classes[][200],
const int* events, const double* weights)
{
char evtitle[200], wttitle[200];
strcpy(evtitle,title);
strcpy(wttitle,title);
strcat(evtitle,": Events");
strcat(wttitle,": Weights");
TCanvas *c
= new TCanvas(canvas,"SPR Input Classes",200,10,600,400);
gStyle->SetPalette(1);
int maxEv = TMath::MaxElement(nClasses,events);
double maxWt = TMath::MaxElement(nClasses,weights);
TPad* pad1 = new TPad("events", evtitle,0,0,1.,0.5);
TPad* pad2 = new TPad("weights",wttitle,0,0.5,1.,1.);
pad1->Draw();
pad2->Draw();
// events
pad1->cd();
TH1I* hev = new TH1I("events",evtitle,nClasses,0,nClasses);
for( int i=0;i<nClasses;i++ )
hev->Fill(classes[i],events[i]);
hev->LabelsDeflate("X");
hev->SetLabelSize(0.06,"X");
hev->SetLabelSize(0.1,"X");
hev->SetLabelSize(0.1,"Y");
hev->SetLineColor(4);
hev->SetFillColor(4);
hev->SetBarWidth(0.8);
hev->SetBarOffset(0.1);
TAxis* yaxis1 = hev->GetYaxis();
yaxis1->SetRangeUser(0.,1.1*maxEv);
hev->Draw("B");
// weights
pad2->cd();
TH1D* hwt = new TH1D("weights",wttitle,nClasses,0,nClasses);
for( int i=0;i<nClasses;i++ )
hwt->Fill(classes[i],weights[i]);
hwt->LabelsDeflate("X");
hwt->SetLabelSize(0.06,"X");
hwt->SetLabelSize(0.1,"X");
hwt->SetLabelSize(0.1,"Y");
hwt->SetLineColor(3);
hwt->SetFillColor(3);
hwt->SetBarWidth(0.8);
hwt->SetBarOffset(0.1);
TAxis* yaxis2 = hwt->GetYaxis();
yaxis2->SetRangeUser(0.,1.1*maxWt);
hwt->Draw("B");
}
示例3: TCanvas
void A1_Mod::PlotCoinc( unsigned aChannelFirst, unsigned aChannelSecond ){
if (fChain == 0) return;
if( c1 ) delete c1;
c1 = new TCanvas( 3 );
stringstream stmp;
unsigned aModule = 0;
unsigned ChannelCoinc = 30;
stmp << "Coincidence - Channels " << aChannelFirst << " : " << aChannelSecond << ";Time [0.1 ns]; Counts";
TH1I * hChannel = new TH1I( "hChannel", stmp.str().c_str(), 10000, 0, 10000 );
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry = 0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry);
nbytes += nb;
// if (Cut(ientry) < 0) continue;
vector<unsigned int> * HITS_PER_CHANNEL = HitsPerChannel[ aModule ];
vector<unsigned int> * COINC = _Data[ aModule ][ ChannelCoinc ];
vector<unsigned int> * TRIG = _Data[ aModule ][ 31 ];
//std::cout << "HITS_PER_CHANNEL 32: " << HITS_PER_CHANNEL->at( 32 ) << std::endl;
if( (HITS_PER_CHANNEL->at( aChannelFirst ) > 0 ) && (HITS_PER_CHANNEL->at( aChannelSecond ) > 0 ) && (HITS_PER_CHANNEL->at( ChannelCoinc ) > 0 ) ){
long ch = COINC->at( 0 );
long trig = TRIG->at( 0 );
std::cout << "Trig - = " << trig - ch << std::endl;
hChannel->Fill( trig - ch );
}
}
hChannel->Draw();
}
示例4: NewCosmicstest
void NewCosmicstest(){
//gROOT->Reset();
TStopwatch *clock0 = new TStopwatch();
TFile *fout = new TFile("Cosmictest.root","RECREATE");
bool high = 1;
//Float_t z0 = 1400;
Float_t z0 = 0;
Float_t yTop = 600;
Float_t xTop = 300;
Float_t zTop = 3650;
Float_t xdist = 3000;
Float_t zdist = 9000;
TH2D *StartXZ = new TH2D("xz","xz;x[cm];z[cm]",30,-(xdist/2),xdist/2,90,z0 - 4500,z0 + 4500);
TH1D *StartTheta = new TH1D("#theta","#theta; #theta_Zenith",100,0,2);
TH1D *StartPhi = new TH1D("#phi","#phi; #phi",50,0,7);
TH1D *StartPHigh = new TH1D("PHigh","P; P[GeV]",100,-1,3);
TH1D *StartP = new TH1D("P","P;P[GeV]",100,-1,3);
TH1D *StartPLow = new TH1D("PLow","P;P[GeV]",100,-1,3);
//TH1D *StartP = new TH1D("P","P;P[GeV]",100,0.1,1000);
BinLogX(StartP);
BinLogX(StartPHigh);
BinLogX(StartPLow);
TH2D *StartPTheta = new TH2D("P,Theta","P-Theta;P[GeV];#theta",150,-1,3,50,0,2);
BinLogX(StartPTheta);
TH2D *MCXZ = new TH2D("MCxz","xz;x[cm];z[cm]",30,-(xdist/2),xdist/2,90,z0 - 4500,z0 + 4500);
TH1D *MCTheta = new TH1D("MC#theta","#theta; #theta",100,0,2);
TH1D *MCPhi = new TH1D("MC#phi","#phi, #phi",50,0,7);
TH1I *MCnTry = new TH1I("nTry","nTry",100,4.6,5);
Double_t totalweightsum=0;
Double_t px,py,pz,x,y,z,w, weighttest, weight;
Int_t nTry,nInside,nEvent,nTest;
Co3Rng *fRandomEngine = new Co3Rng();
int EVENTS = 400000;
Int_t kmax = 40000;
float weight1,weight1Low,weight1High;
weight1Low = 123*xdist*zdist/EVENTS/10000; // expected #muons per spill/ #simulated events per spill 174*30*90/500000
cout<<weight1Low<<endl;
double I = fRandomEngine->fSpectrumH->Integral(100,1000);
weight1High = 2*TMath::Pi()/3*I*xdist*zdist/EVENTS/10000;
//weight2 = 900/I; // 1/(mean momentum weight), P_max-P_min/(3*0.3044/2pi)
cout<< weight1High<<endl;
Float_t weight3 = 4.833931503; // MC average of nTry/nEvents 4.833949997 +- 0.000010494
weight1 = 1;
weight = weight1 / weight3;
nInside = 0; nEvent = 0; nTest = 0; weighttest = 0;
y = 1900; //20m over beam axis
w = weight/kmax;
clock0->Start();
for(Int_t k = 0;k<kmax;k++){
cout<<k<<endl;
nTry =0;
for(Int_t i=0;i<EVENTS;i++){
Bool_t hit = 0;
do{
// shower characteristics
double phi = fRandomEngine->Uniform(0,2*TMath::Pi());
double theta = fRandomEngine->fTheta->GetRandom();
//momentum components
px = TMath::Sin(phi)*TMath::Sin(theta);
pz = TMath::Cos(phi)*TMath::Sin(theta);
py = -TMath::Cos(theta);
// start position, area 1120m^2
x = fRandomEngine->Uniform(-xdist/2,xdist/2);
z = fRandomEngine->Uniform(z0 - zdist/2, z0 + zdist/2);
// claim for flight close to the actual detector
if((abs(x-(y+yTop)*px/py) < xTop && abs(z-z0-(y+yTop)*pz/py) < zTop) || (abs(x-(y-yTop)*px/py) < xTop && abs(z-z0-(y-yTop)*pz/py) < zTop)|| abs(y-(x+xTop)*py/px)<yTop && abs(z-z0-(x+xTop)*pz/px)<zTop || abs(y-(x-xTop)*py/px)<yTop && abs(z-z0-(x-xTop)*pz/px)<zTop){
// muon momentum
double P;
//if (!high) {P = fRandomEngine->NEWstest(theta);}
//else {P = fRandomEngine->fSpectrumH->GetRandom();}
// high
P = fRandomEngine->fSpectrumH->GetRandom();
w = weight1High/weight3/kmax;
StartP->Fill(P,w);
StartPHigh->Fill(P,w);
StartPTheta->Fill(P,theta,w);
// low
P = fRandomEngine->NEWstest(theta);
w = weight1Low/weight3/kmax;
StartP->Fill(P,w);
StartPLow->Fill(P,w);
StartPTheta->Fill(P,theta,w);
//.........这里部分代码省略.........
示例5: main
int main(int argc, char* argv[])
{
TFile* file = new TFile("neutrons.root", "recreate");
TChain* chain = createChain(argc, argv);
// Assign addresses
double target_charge, veto_charge, target_cb, veto_cb;
unsigned long long time;
chain->SetBranchAddress("target_total", &target_charge);
chain->SetBranchAddress("veto_total", &veto_charge);
chain->SetBranchAddress("target_cb", &target_cb);
chain->SetBranchAddress("veto_cb", &veto_cb);
chain->SetBranchAddress("time", &time);
const int chainEntries = chain->GetEntries();
const double adc2us = 4/1000.0;
const int maxCharge = 30000;
TH1F* histogram = new TH1F("timing", "timing", 2000, 0, 2000);
TH1F* neutrons = new TH1F("nspec", "neutron spectrum", 200, 0, maxCharge);
TH1F* hydrogen = new TH1F("nhspec", "capture on hydrogen", 200, 0, maxCharge);
TH1F* externals = new TH1F("externspec", "externals", 200, 0, maxCharge);
TH1F* tsa_neutrons = new TH1F("tsa_spec", "neutron spectrum from tsa plot", 200, 0, maxCharge);
TH1I* multiplicity = new TH1I("mult", "Event multiplicity", 20, 0, 20);
unsigned long long t1=0;
unsigned long long t2=0;
unsigned long long t3=0;
// cuts
const double max_target_cb = 0.5;
const double max_veto_charge = 500;
const double min_target_charge = 0;
const double tale_min = 800;
const double tale_max = 2000;
const double neutrons_min = 0;
const double neutrons_time = 50;
const double hydro_min = 100;
const double hydro_max = 150;
double previous_charge = 0;
// TSA Plots
const int nbins = 100;
double logmin = -0.3;
double logmax = 4.5;
double binwidth = (logmax-logmin)/double(nbins);
double bin_array[nbins+1];
bin_array[0]=std::pow(10, logmin);
for(int i=1; i<=nbins; i++)
bin_array[i]=bin_array[0]+std::pow(10,logmin + i*binwidth);
TH2F* tsa = new TH2F("tsa", "time series analysis", nbins, bin_array, nbins, bin_array);
// Keep track of timing in order to get multiplicity
std::vector<unsigned long long> time_vec;
time_vec.resize(1, 0);
for(int evt=0; evt<chainEntries; ++evt)
{
// Print out progress
if(!(evt%1000))
std::cout << std::floor(double(evt)/chainEntries*100) << "%\r";
previous_charge = target_charge;
chain->GetEvent(evt);
if(veto_charge<max_veto_charge &&
target_charge>min_target_charge &&
target_cb < max_target_cb )
{
t3=t2;
t2=t1;
t1=time;
double d12 = (t1-t2)*adc2us;
double d23 = (t2-t3)*adc2us;
tsa->Fill(d23, d12);
if((time-time_vec[0])*adc2us < 50)
time_vec.push_back(time);
else
{
multiplicity->Fill(time_vec.size());
time_vec.clear();
time_vec.push_back(time);
}
histogram->Fill(d12);
if(d12 < neutrons_time && d12 > neutrons_min)
neutrons->Fill(target_charge);
if(d12 < tale_max && d12 > tale_min)
externals->Fill(target_charge);
if(d12 < hydro_max && d12 > hydro_min)
hydrogen->Fill(target_charge);
if(d12 < neutrons_time && d23 < neutrons_time)
tsa_neutrons->Fill(previous_charge);
}
}
std::cout << std::endl;
// Events in tale:
TF1* fitter_extern = new TF1("fitter_extern", "[0]*TMath::Exp([1]*x)",
0, 2000);
fitter_extern->SetParameters(1000, -4e-4);
//.........这里部分代码省略.........
示例6: processWaveforms
void processWaveforms(TTree* treeToSort, vector<Plots>& targetPlots, TFile* outFile, string mode)
{
TH1D* fittedTimeHisto = new TH1D("raw fitted times", "raw fitted times", TOF_BINS,TOF_LOWER_BOUND,TOF_RANGE);
TH2D* deltaTVsPulseHeight = new TH2D("delta T vs. pulse height","delta T vs. pulse height", 200, -10, 10, 1580, 0, 15800);
TH2D* deltaTVsPulseIntegral0 = new TH2D("delta T vs. pulse integral, target 0","delta T vs. pulse integral, target 0",100, -10, 10, pow(2,15), 0, pow(2,15));
TH2D* deltaTVsPulseIntegral1 = new TH2D("delta T vs. pulse integral, target 1","delta T vs. pulse integral, target 1",100, -10, 10, pow(2,15), 0, pow(2,15));
TH2D* deltaTVsPulseIntegral2 = new TH2D("delta T vs. pulse integral, target 2","delta T vs. pulse integral, target 2",100, -10, 10, pow(2,15), 0, pow(2,15));
TH2D* deltaTVsPulseIntegral3 = new TH2D("delta T vs. pulse integral, target 3","delta T vs. pulse integral, target 3",100, -10, 10, pow(2,15), 0, pow(2,15));
TH2D* deltaTVsPulseIntegral4 = new TH2D("delta T vs. pulse integral, target 4","delta T vs. pulse integral, target 4",100, -10, 10, pow(2,15), 0, pow(2,15));
TH2D* deltaTVsPulseIntegral5 = new TH2D("delta T vs. pulse integral, target 5","delta T vs. pulse integral, target 5",100, -10, 10, pow(2,15), 0, pow(2,15));
TH1D* triggerAmplitudeHisto = new TH1D("triggerAmplitudeHisto","triggerAmplitudeHisto",pow(2,14),0,pow(2,14));
TH1D* relativeTriggerTimeHisto = new TH1D("relativeTriggerTimeHisto","relative trigger time, from start of fitted wavelet",200,-5,5);
TH2D* relativeTriggerTimeVsAmplitude = new TH2D("relativeTriggerTimeVSAmplitude","relative trigger time vs amplitude", 200, -5, 5, pow(2,14), 0, pow(2,14));
TH1D* gammaToGammaTimeH = new TH1D("gammaToGammaTimeH","time between consecutive gammas", 1000, -5, 5);
if(mode=="DPP")
{
setBranchesHistos(treeToSort);
int totalEntries = treeToSort->GetEntries();
cout << "Total waveforms = " << totalEntries << endl;
vector<double> triggerList;
waveformWrap = new TMultiGraph("DPP waveforms", "DPP waveforms");
vector<TGraph*> waveletGraphs;
vector<TGraph*> triggerGraphs;
int gammaGate[2] = {80,90};
double prevGammaTime = 0;
for(int j=1; j<totalEntries; j++)
{
if(j%1000==0)
{
cout << "Processing triggers on waveform " << j << "\r";
fflush(stdout);
}
/*if(j>500)
{
break;
}*/
triggerList.clear();
triggerValues.clear();
// pull individual waveform event
treeToSort->GetEntry(j);
// calculate the baseline for this waveform
BASELINE = calculateBaseline(*procEvent.waveform);
// Loop through all points in the waveform and fit peaks
for(int k=DPP_PEAKFIT_START; (size_t)k<procEvent.waveform->size(); k++)
{
// Check to see if this point creates a new trigger
if(isTrigger(k, *procEvent.waveform))
{
// trigger found - plot/fit/extract time
double timeDiff = procEvent.completeTime-procEvent.macroTime;
double microTime = fmod(timeDiff,MICRO_LENGTH);
if(microTime > gammaGate[0] && microTime < gammaGate[1])
{
processTrigger(j, k, triggerList, *procEvent.waveform);
double fullTime = procEvent.completeTime-procEvent.macroTime+data.trigger1Time+DPP_PEAKFIT_OFFSET;
gammaToGammaTimeH->Fill(fmod(fullTime,MICRO_LENGTH)-prevGammaTime);
prevGammaTime=fmod(fullTime,MICRO_LENGTH);
fillTriggerHistos(fullTime, targetPlots);
fittedTimeHisto->Fill(data.trigger1Time);
}
//processTrigger(j, k, triggerList, *procEvent.waveform);
break;
}
}
//produceTriggerOverlay(j, triggerList, *procEvent.waveform);
TH2D* deltaTVsPulseIntegralHisto;
switch(procEvent.targetPos-1)
{
case 0:
deltaTVsPulseIntegralHisto = deltaTVsPulseIntegral0;
break;
case 1:
deltaTVsPulseIntegralHisto = deltaTVsPulseIntegral1;
break;
//.........这里部分代码省略.........
示例7: ExtractTrackBasedTiming
//.........这里部分代码省略.........
outFile.close(); // clear file
int nBinsX = thisHist->GetNbinsX();
int nBinsY = thisHist->GetNbinsY();
TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
for (int i = 1 ; i <= nBinsX; i++){
TH1D *projY = thisHist->ProjectionY("temp", i, i);
// Scan over the histogram
//chose the correct number of bins based on the histogram
float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
float timeWindow = 3; //ns (Full Width)
int binWindow = int(timeWindow / nsPerBin);
double maxEntries = 0;
double maxMean = 0;
for (int j = 1 ; j <= projY->GetNbinsX();j++){
int minBin = j;
int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
double sum = 0, nEntries = 0;
for (int bin = minBin; bin <= maxBin; bin++){
sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
nEntries += projY->GetBinContent(bin);
if (bin == maxBin){
if (nEntries > maxEntries) {
maxMean = sum / nEntries;
maxEntries = nEntries;
}
}
}
}
//In the case there is RF, our job is to pick just the number of the correct beam bunch, so that's really all we need.
if(useRF) {
int beamBucket = int((maxMean / RF_Period) + 0.5); // +0.5 to handle rounding correctly
selectedTAGMOffset->SetBinContent(i, beamBucket);
TAGMOffsetDistribution->Fill(beamBucket);
}
else{
selectedTAGMOffset->SetBinContent(i, maxMean);
TAGMOffsetDistribution->Fill(maxMean);
}
}
double meanOffset = TAGMOffsetDistribution->GetMean();
// This might be in units of beam bunches, so we need to convert
if (useRF) meanOffset *= RF_Period;
if (verbose) {
cout << "Dumping TAGM results...\n=======================================" << endl;
cout << "TAGM mean Offset = " << meanOffset << endl;
cout << "fADC Offsets" << endl;
}
outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out);
//for (int i = 1 ; i <= nBinsX; i++){
// Loop over rows
if (verbose) cout << "Column\tRow\tvalueToUse\toldValue\tmeanOffset\tTotal" << endl;
for (unsigned int column = 1; column <= 102; column++){
int index = GetCCDBIndexTAGM(column, 0);
double valueToUse = selectedTAGMOffset->GetBinContent(index);
if (useRF) valueToUse *= RF_Period;
//if (valueToUse == 0) valueToUse = meanOffset;
outFile << "0 " << column << " " << valueToUse + tagm_fadc_time_offsets[index-1] - meanOffset<< endl;
if (verbose) printf("0\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", column, valueToUse, tagm_fadc_time_offsets[index-1], meanOffset,
valueToUse + tagm_fadc_time_offsets[index-1] - meanOffset);
if (column == 9 || column == 27 || column == 81 || column == 99){
for (unsigned int row = 1; row <= 5; row++){
index = GetCCDBIndexTAGM(column, row);
valueToUse = selectedTAGMOffset->GetBinContent(index);
示例8: runcorr
void runcorr(int condor_iter, int trackqual, string flist = "", string tag = "", int nmin = 220, int nmax = 1000, float pttrigmin = 1, float pttrigmax = 2, float ptassmin = 1, float ptassmax = 1)
{
const int nptbins = 1;
const int ncentbins = 1;
const int najbins = 1;
string buffer;
vector<string> listoffiles;
int nlines = 0;
// ifstream infile("sortedforests.txt");
// ifstream infile("doeproposalforests.txt");
// ifstream infile("PA2013_PromptReco_Json_FullTrack12_v84_prod2.txt");
// ifstream infile("PA2013_PromptReco_Json_FullTrack12_v84_prod2_MIT.txt");
// ifstream infile("HIMinBiasUPC_PbPbUPC_pptracking_452p1_forest_v2.txt");
ifstream infile(flist.data());
if (!infile.is_open()) {
cout << "Error opening file. Exiting." << endl;
return;
} else {
while (!infile.eof()) {
infile >> buffer;
listoffiles.push_back(buffer);
nlines++;
}
}
bool dostdhists = (condor_iter%(nptbins * ncentbins * najbins) == 0);
int ptbin = condor_iter % nptbins;
int centbin = (condor_iter / nptbins) % ncentbins;
int ajbin = (condor_iter / (nptbins * ncentbins)) % najbins;
int filenum = (condor_iter / (nptbins * ncentbins * najbins));
cout << "ipt: " << ptbin << " icent: " << centbin << " iaj: " << ajbin << " filenum: " << filenum << " dostdhists: " << dostdhists << endl;
// int nmin = 220 , nmax = 1000, one = 1;
//! for first iteration of forest production
/*
if(filenum==0) { nmin = 110 ; nmax = 1000; }
if(filenum==1) { nmin = 90 ; nmax = 110 ; }
if(filenum>1 ) { nmin = 0 ; nmax = 35 ; }
if(filenum>9 ) { nmin = 35 ; nmax = 90 ; }
*/
//! for second iteration of forest production
// /*
// if(filenum<24 ) { nmin = 90 ; nmax = 110 ; }
// if(filenum<22 ) { nmin = 35 ; nmax = 90 ; }
// if(filenum<12 ) { nmin = 110 ; nmax = 1000; }
// if(filenum<10 ) { nmin = 0 ; nmax = 35 ; }
// */
//! for second iteration of forest production
/*
if(filenum<26 ) { nmin = 90 ; nmax = 110 ; }
if(filenum<23 ) { nmin = 35 ; nmax = 90 ; }
if(filenum<13) { nmin = 110 ; nmax = 1000; }
if(filenum<10) { nmin = 0 ; nmax = 35 ; }
*/
corrana(listoffiles[filenum].data(),trackqual);
// double pttriglow[] = {1 ,14 ,16 ,20, 1};
// double pttrighigh[] = {2 ,16 ,20 ,24, 3};
// double ptasslow[] = {1 ,1 ,1 ,1 , 1};
// double ptasshigh[] = {2 ,2 ,2 ,2 , 3};
int centmin[] = {0,4,8,12,16,20,24,28,32};
int centmax[] = {41,8,12,16,20,24,28,32,36};
// TFile * outf = new TFile(Form("%s_trkqaul%d_nmin%d_nmax%d_tptmin%d_tptmax%d_aptmin%d_aptmax%d_%d.root",tag.data(),trackqual,nmin,nmax,(int)pttrigmin/one,(int)pttrigmax/one,(int)ptassmin/one,(int)ptassmax/one,filenum),"recreate");
cout<<Form("%s_trkqaul%d_nmin%d_nmax%d_tptmin%d_tptmax%d_aptmin%d_aptmax%d_%d.root",tag.data(),trackqual,nmin,nmax,(int)pttrigmin,(int)pttrigmax,(int)ptassmin,(int)ptassmax,filenum)<<endl;
TFile * outf = new TFile(Form("%s_trkqaul%d_nmin%d_nmax%d_tptmin%d_tptmax%d_aptmin%d_aptmax%d_%d.root",tag.data(),trackqual,nmin,nmax,(int)pttrigmin,(int)pttrigmax,(int)ptassmin,(int)ptassmax,filenum),"recreate");
// int i = 0;
int cent = 0;
cout<<"cent iteration "<<cent<<endl;
TH2D * ttsig = TrackTrackSignal(pttrigmin,pttrigmax,ptassmin,ptassmax,centmin[cent],centmax[cent],nmin,nmax);
TH2D * ttbak = TrackTrackBackground(pttrigmin,pttrigmax,ptassmin,ptassmax,centmin[cent],centmax[cent],nmin,nmax);
TH1I * hntottrig = new TH1I(Form("nttottrig_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttrigmin,(int)pttrigmax,(int)ptassmin,(int)ptassmax,centmin[cent],centmax[cent]),"",1,0.5,1.5);
int myntottrig = GetNTotTrig();
hntottrig->Fill(1,myntottrig);
cout<<"ntottrig: "<<myntottrig<<endl;
TH2D * ttcorr = (TH2D*)ttsig->Clone(Form("corr_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttrigmin,(int)pttrigmax,(int)ptassmin,(int)ptassmax,centmin[cent],centmax[cent]));
ttcorr->Divide(ttbak);
ttcorr->Scale(ttbak->GetBinContent(ttbak->FindBin(0,0)));
ttcorr->Scale(1/0.0594998609); //! bin width
ttcorr->GetXaxis()->SetRange(ttcorr->GetXaxis()->FindBin(-4.0),ttcorr->GetXaxis()->FindBin(4.0));
ttcorr->GetYaxis()->SetRange(ttcorr->GetYaxis()->FindBin(-3.1415926/2.0),ttcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
outf->Write();
outf->Close();
}
示例9: ExtractTrackBasedTiming
//.........这里部分代码省略.........
outFile.close(); // clear file
int nBinsX = thisHist->GetNbinsX();
int nBinsY = thisHist->GetNbinsY();
TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
for (int i = 1 ; i <= nBinsX; i++){
TH1D *projY = thisHist->ProjectionY("temp", i, i);
// Scan over the histogram
//chose the correct number of bins based on the histogram
float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
float timeWindow = 3; //ns (Full Width)
int binWindow = int(timeWindow / nsPerBin);
double maxEntries = 0;
double maxMean = 0;
for (int j = 1 ; j <= projY->GetNbinsX();j++){
int minBin = j;
int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
double sum = 0, nEntries = 0;
for (int bin = minBin; bin <= maxBin; bin++){
sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
nEntries += projY->GetBinContent(bin);
if (bin == maxBin){
if (nEntries > maxEntries) {
maxMean = sum / nEntries;
maxEntries = nEntries;
}
}
}
}
//In the case there is RF, our job is to pick just the number of the correct beam bunch, so that's really all we need.
if(useRF) {
int beamBucket = int((maxMean / RF_Period) + 0.5); // +0.5 to handle rounding correctly
selectedTAGMOffset->SetBinContent(i, beamBucket);
TAGMOffsetDistribution->Fill(beamBucket);
}
else{
selectedTAGMOffset->SetBinContent(i, maxMean);
TAGMOffsetDistribution->Fill(maxMean);
}
}
/*
if (!useRF){
//TFitResultPtr fr1 = selectedTAGMOffset->Fit("pol1", "SQ", "", 0.5, nBinsX + 0.5);
TFitResultPtr fr1 = selectedTAGMOffset->Fit("pol1", "SQ", "", 5, 50);
for (int i = 1 ; i <= nBinsX; i++){
double x0 = fr1->Parameter(0);
double x1 = fr1->Parameter(1);
//double x2 = fr1->Parameter(2);
//double fitResult = x0 + i*x1 + i*i*x2;
double fitResult = x0 + i*x1;
double outlierCut = 20;
double valueToUse = selectedTAGMOffset->GetBinContent(i);
if (fabs(selectedTAGMOffset->GetBinContent(i) - fitResult) > outlierCut && valueToUse != 0.0){
valueToUse = fitResult;
}
selectedTAGMOffset->SetBinContent(i, valueToUse);
if (valueToUse != 0 ) TAGMOffsetDistribution->Fill(valueToUse);
}
}
*/
double meanOffset = TAGMOffsetDistribution->GetMean();
// This might be in units of beam bunches, so we need to convert
if (useRF) meanOffset *= RF_Period;
示例10: ExtractTimeOffsetsAndCEff
//.........这里部分代码省略.........
double globalOffset[768];
TH1D * selectedBCALOffset = new TH1D("selectedBCALOffset", "Selected Global BCAL Offset; CCDB Index; Offset [ns]", 768, 0.5, 768 + 0.5);
TH1I * BCALOffsetDistribution = new TH1I("BCALOffsetDistribution", "Global BCAL Offset; Global Offset [ns]; Entries", 100, -10, 10);
thisHist = Get2DHistogram("BCAL_Global_Offsets", "Target Time", "Target Time Minus RF Time Vs. Cell Number");
if(thisHist != NULL) {
int nBinsX = thisHist->GetNbinsX();
int nBinsY = thisHist->GetNbinsY();
for (int i = 1 ; i <= nBinsX; i++) {
TH1D *projY = thisHist->ProjectionY("temp", i, i);
// Scan over the histogram
float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
float timeWindow = 0.5; //ns (Full Width)
int binWindow = int(timeWindow / nsPerBin);
double maxEntries = 0;
double maxMean = 0;
for (int j = 1 ; j <= projY->GetNbinsX(); j++) {
int minBin = j;
int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
double sum = 0, nEntries = 0;
for (int bin = minBin; bin <= maxBin; bin++) {
sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
nEntries += projY->GetBinContent(bin);
if (bin == maxBin) {
if (nEntries > maxEntries) {
maxMean = sum / nEntries;
maxEntries = nEntries;
}
}
}
}
globalOffset[i-1] = maxMean;
selectedBCALOffset->SetBinContent(i, maxMean);
BCALOffsetDistribution->Fill(maxMean);
}
}
outputFile->cd("Fits");
// Now we want to loop through all available module/layer/sector and try to make a fit of each one
for (unsigned int iModule = 1; iModule <=48; iModule++) {
for (unsigned int iLayer = 1; iLayer <= 4; iLayer++) { // Only 3 layers with TDCs
for (unsigned int iSector = 1; iSector <= 4; iSector++) {
int the_cell = (iModule - 1) * 16 + (iLayer - 1) * 4 + iSector;
int the_tdc_cell = (iModule - 1) * 12 + (iLayer - 1) * 4 + iSector; // One less layer of TDCs
// Format the string to lookup the histogram by name
char name[200];
sprintf(name, "Module %.2i Layer %.2i Sector %.2i", iModule, iLayer, iSector);
// These histograms are created on the fly in the plugin, so there is a chance that they do not exist, in which case the pointer will be NULL
TH2I *h_offsets = Get2DHistogram ("BCAL_TDC_Offsets", "Z Position", name);
// Use FitSlicesY routine to extract the mean of each x bin
TObjArray ySlices;
if (h_offsets != NULL) {
h_offsets->RebinX(5);
TProfile *profile = h_offsets->ProfileX();
f1->SetParameters(0, 1); // Just out initial guess
TFitResultPtr fr = profile->Fit(f1, "SQR");
Int_t fitStatus = fr;
if (fitStatus == 0) {
double c0 = fr->Parameter(0);
double c0_err = fr->ParError(0);
double c1 = fr->Parameter(1);
double c1_err = fr->ParError(1);
示例11: Loop
void Zlumi::Loop()
{
// In a ROOT session, you can do:
// Root > .L Zlumi.C
// Root > Zlumi t
// Root > t.GetEntry(12); // Fill t data members with entry number 12
// Root > t.Show(); // Show values of entry 12
// Root > t.Show(16); // Read and show values of entry 16
// Root > t.Loop(); // Loop on all entries
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus("branchname",1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
gROOT->ForceStyle();
tdrStyle();
if (fChain == 0) return;
int minRun=0;
int maxRun=0;
int maxLS=0;
bool forminRun=true;
//TH2I * LumiSRun = new TH2I("LumiSRun", "LS vs Run", 3000, 0., 3000., 20000, 160000., 180000.);
//TH1I * test2 = new TH1I("test2","test2", 3000,0, 3000);
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
if(forminRun && (Run!=0)){minRun=Run; forminRun=false;}
if((Run!=0) && (Run>maxRun)){maxRun=Run;}
if((LS!=0) && (LS>maxLS)){maxLS=LS;}
//printf("run %i ls %i \n",Run,LS);
}
cout << nentries << " nentries \n";
TH2I *LumiSRun = new TH2I("LumiSRun", "LS vs Run", maxLS, 0, maxLS, maxRun-minRun+2, minRun-1, maxRun+1);
TH1I *Runs = new TH1I("Runs","Run", maxRun-minRun+2, minRun-1, maxRun+1);
Runs->Sumw2();
nbytes = 0; nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
LumiSRun->Fill(LS,Run);
//printf("run %i ls %i \n",Run,LS);
Runs->Fill(Run);
//test2->Fill(LS);
}
printf("minRun %i maxRun %i \n",minRun,maxRun);
LumiSRun->Draw();
for (int h=0;h<Runs->GetNbinsX();h++){
Runs->SetBinError(h+1,sqrt(Runs->GetBinContent(h+1)) );
}
Runs->Draw();
//test->Draw();
//test2->Draw();
TH1F *FileRuns = new TH1F("FileRuns","Run from Lumicalc", maxRun-minRun+2, minRun-1, maxRun+1);
TH1D *XsecDistro = new TH1D("XsecDistro","X sec distribution", 60, 0., 0.6);
//-------------
const Int_t mpt = maxRun-minRun;
int fileRun[mpt];
double Lumi[mpt];
int npt = 0;
// read data file
ifstream file;
//file.open("./2011-run-lumi.txt");
file.open("./LumiAeB-dav.txt");
while (1) {
file >> fileRun[npt] >> Lumi[npt];
if ( ! file.good() ) break;
cout << "x = " << fileRun[npt] << " y = " << Lumi[npt] << endl;
FileRuns->SetBinContent((fileRun[npt]-minRun+2),0.2);
npt++;
//.........这里部分代码省略.........