Skip to content

Instantly share code, notes, and snippets.

@jyoshida-sci
Last active January 18, 2022 21:42
Show Gist options
  • Save jyoshida-sci/9a1f69fb3fcdda5d6cd6 to your computer and use it in GitHub Desktop.
Save jyoshida-sci/9a1f69fb3fcdda5d6cd6 to your computer and use it in GitHub Desktop.
ROOT macros. validated by ROOT 5.34/24

#ROOT Macros

##TNtuple_TH1D.c How to make a 2d-scatter-plot and a 1d-projected-histogram.

##pyroot_demo.py How to draw a graph of f(x) = abs(sin(x)/x).

##randomnumber.c Random number generator.

##pol1fit.c linear fitting for pol1.txt; y=3x+1.

##pol6fit.c polynomial curve fitting for pol6.txt made with makepol6.c.

##rootmacro_args.c How to pass args from commandline.

##drawArrows.c How to draw arrows on TH2D.

##fillToTH1D.c How to fill and merge data into a TH1D with TNtuple::Draw().

##TH1D_addnormalize.c Filling gaussian-random-number in a TH1D, merging 2 TH1Ds and normalizing.

##TNtuple_TCut.c Alias and TCut for TNtuple.

##changing_cutcond.c How to change a cut condition with sprintf and TCut. data.txt by gen_data.py

##derivative2.c and derivative2.py Inflection points of tanh(f*gaus(x,mean,sigma)) as x where second-derivative=0.

##chi2distribution.c How to generate Chi-squared distribution (n.d.f=3).

##chi2table.c How to generate a table of upper critical values of chi-square distribution.

##poissonhist.c Histogram of Poisson distribution.

##TDatabasePDG.c How to use TDatabasePDG.

##rms_mean.c TMath::Mean and TMath::RMS.

##meanFreePathOfKminus.c Mean free path of Kminus meson in nuclear emulsion as a function of laboratory energy

##decayEventGen.c

1 11
2 22
3 33
4 44
//usage: >root changing_cutcond.c
{
TCanvas* c1 = new TCanvas("c1","canvas",800,800);
TNtuple* nt = new TNtuple("nt","data","i:j:r");
nt->SetMarkerStyle(8);
nt->ReadFile("data.txt");
nt->Draw("i:j:r","","colz");
c1->Print("a.png");
for(int i=1; i<7; i++){
char condition[128];
sprintf(condition,"r<%d",i);
printf("%d %s\n",i,condition);
TCut mycut = condition;
nt->Draw("i:j:r",mycut,"colz");
char filename[128];
sprintf(filename,"i_j_r%02d.png",i);
c1->Print(filename);
}
}
//usage: >root chi2distribution.c
Double_t myfunc(Double_t *x,Double_t *par)
{
Float_t xx =x[0];
Double_t f = ROOT::Math::chisquared_pdf(xx, par[0]);
return f;
}
chi2distribution(){
const int NDF=3;
const int ROUND=1E5;
TH1D* h = new TH1D("h","Chi-squared distribution",14, 0.0, 14.0);
for(int n=0; n<ROUND; n++){
double chisq = 0.0;
for(int i=0; i<NDF; i++){
double sigma = gRandom->Uniform(0.5, 2.0);
double r = gRandom->Gaus(0.0, sigma);
chisq += (r*r) / (sigma*sigma);
}
h->Fill(chisq);
}
//normalize
h->Scale( 1.0 / h->GetEntries());
h->SetMaximum(0.3);
h->SetMinimum(0.0);
h->Draw();
//range(0-40.0) number_of_parameter:1
TF1* pdf = new TF1("pdf", myfunc, 0.0, 40.0, 1);
pdf->SetParameter(0,NDF);
pdf->Draw("same");
}
//usage >root this.c
{
double prob[12] = {0.995, 0.990, 0.975, 0.950, 0.050, 0.025, 0.010, 0.005};
printf("Upper critical values of chi-square distribution\n");
for(int i=0; i<8; i++){
printf("%.4f ",prob[i]);
}
printf("\n==============================================================\n");
for(int ndf=1; ndf<20; ndf++){
for(int i=0; i<8; i++){
double val = ROOT::Math::chisquared_quantile_c(prob[i], ndf);
printf("%.4f ",val);
}
printf("\n");
}
cout << "======" << endl;
cout << ROOT::Math::chisquared_quantile_c(0.995, 5) << endl;
cout << ROOT::Math::chisquared_quantile(1.0-0.995, 5) << endl;
cout << TMath::Prob(0.411742, 5) << endl;
cout << ROOT::Math::chisquared_cdf_c(0.411742, 5) << endl;
cout << 1.0 - ROOT::Math::chisquared_cdf(0.411742, 5) << endl;
}
/*
Upper critical values of chi-square distribution
0.9950 0.9900 0.9750 0.9500 0.0500 0.0250 0.0100 0.0050
==============================================================
0.0000 0.0002 0.0010 0.0039 3.8415 5.0239 6.6349 7.8794
0.0100 0.0201 0.0506 0.1026 5.9915 7.3778 9.2103 10.5966
0.0717 0.1148 0.2158 0.3518 7.8147 9.3484 11.3449 12.8382
0.2070 0.2971 0.4844 0.7107 9.4877 11.1433 13.2767 14.8603
0.4117 0.5543 0.8312 1.1455 11.0705 12.8325 15.0863 16.7496
0.6757 0.8721 1.2373 1.6354 12.5916 14.4494 16.8119 18.5476
0.9893 1.2390 1.6899 2.1673 14.0671 16.0128 18.4753 20.2777
1.3444 1.6465 2.1797 2.7326 15.5073 17.5345 20.0902 21.9550
1.7349 2.0879 2.7004 3.3251 16.9190 19.0228 21.6660 23.5894
2.1559 2.5582 3.2470 3.9403 18.3070 20.4832 23.2093 25.1882
2.6032 3.0535 3.8157 4.5748 19.6751 21.9200 24.7250 26.7568
3.0738 3.5706 4.4038 5.2260 21.0261 23.3367 26.2170 28.2995
3.5650 4.1069 5.0088 5.8919 22.3620 24.7356 27.6882 29.8195
4.0747 4.6604 5.6287 6.5706 23.6848 26.1189 29.1412 31.3193
4.6009 5.2293 6.2621 7.2609 24.9958 27.4884 30.5779 32.8013
5.1422 5.8122 6.9077 7.9616 26.2962 28.8454 31.9999 34.2672
5.6972 6.4078 7.5642 8.6718 27.5871 30.1910 33.4087 35.7185
6.2648 7.0149 8.2307 9.3905 28.8693 31.5264 34.8053 37.1565
6.8440 7.6327 8.9065 10.1170 30.1435 32.8523 36.1909 38.5823
======
0.411742
0.411742
0.995
0.995
0.995
(int)1922239832
*/
0 0 8.90378653913
0 1 6.71033038022
0 2 6.11916065257
0 3 1.33405828703
0 4 9.59983273829
0 5 9.26842007104
0 6 4.71684647584
0 7 3.90643868036
0 8 0.427981997306
0 9 3.0997863775
0 10 6.82465265906
0 11 0.085964755748
0 12 3.58159648187
0 13 0.919485672947
0 14 6.68038345382
0 15 1.82665711806
0 16 8.84772209217
0 17 0.0131768786963
0 18 4.38540036598
0 19 9.0183199358
1 0 6.26627627465
1 1 4.28616177958
1 2 9.01845780208
1 3 2.28896842619
1 4 2.47511838301
1 5 6.44209071018
1 6 8.72579598619
1 7 7.4695199596
1 8 2.5776714859
1 9 8.19811946262
1 10 7.1022360782
1 11 9.30512805047
1 12 0.980422344685
1 13 5.11269940064
1 14 7.59247528586
1 15 2.0453038446
1 16 8.57134615295
1 17 8.48714367256
1 18 8.89873016598
1 19 7.04013821548
2 0 6.36567052518
2 1 4.63827765187
2 2 1.84021622866
2 3 2.50052995834
2 4 6.01150249332
2 5 7.95536081342
2 6 1.85469916167
2 7 6.22819110157
2 8 8.51771650057
2 9 6.06918713528
2 10 3.59320117284
2 11 0.918993513523
2 12 5.95678596858
2 13 6.19942198613
2 14 0.9440926235
2 15 5.18106912703
2 16 9.22721229073
2 17 0.0418175752421
2 18 0.216802238949
2 19 6.62517238255
3 0 5.19312659232
3 1 0.0793717251225
3 2 8.6419783976
3 3 3.71330983442
3 4 7.59652343361
3 5 6.26304610509
3 6 9.47573411335
3 7 9.85390031243
3 8 4.47108064967
3 9 3.19726461632
3 10 4.68251905522
3 11 9.41324450977
3 12 3.7205538379
3 13 8.2413221967
3 14 1.04407296894
3 15 0.00989152670864
3 16 4.61621539669
3 17 5.5776741532
3 18 4.76933513138
3 19 1.39809527791
4 0 1.80217143173
4 1 7.58771754085
4 2 2.3505179615
4 3 7.17010501795
4 4 2.16368851471
4 5 2.11158474316
4 6 0.132313336339
4 7 2.22912341094
4 8 8.05695819345
4 9 9.70050772229
4 10 9.06763707362
4 11 7.43598828057
4 12 7.65959804334
4 13 0.913972617946
4 14 6.94391160787
4 15 2.10045239613
4 16 6.66169244801
4 17 4.06872387053
4 18 5.97055040238
4 19 1.47555852646
5 0 5.60390805676
5 1 3.25137533817
5 2 0.808574878951
5 3 2.18859886953
5 4 8.35968710666
5 5 5.57236776805
5 6 9.34059953065
5 7 8.97661879376
5 8 5.33716066815
5 9 8.77520502801
5 10 2.54359338582
5 11 6.51334674606
5 12 9.79581562796
5 13 7.89232517279
5 14 9.11985623888
5 15 6.56234063117
5 16 4.72884886441
5 17 2.07115422696
5 18 4.18960861862
5 19 0.142283683922
6 0 8.83781306949
6 1 6.65505236195
6 2 7.27911303565
6 3 6.92496610067
6 4 5.80835140395
6 5 9.57716523268
6 6 9.61031484424
6 7 1.75281646079
6 8 1.91703928637
6 9 8.21967898262
6 10 1.88583406958
6 11 0.464392304938
6 12 1.63616246252
6 13 8.29576779408
6 14 7.26810512715
6 15 5.97583921
6 16 2.02080882697
6 17 2.80760262137
6 18 5.55185501403
6 19 6.07665661864
7 0 7.1092293555
7 1 4.27489152853
7 2 6.72839610797
7 3 9.58114926987
7 4 2.71552640628
7 5 1.67864841929
7 6 4.7502938148
7 7 0.508042510421
7 8 3.13142926694
7 9 7.95496736395
7 10 5.93883209837
7 11 1.74941912841
7 12 0.412688003401
7 13 9.88298493271
7 14 8.32393389942
7 15 4.18482217153
7 16 4.89611164619
7 17 3.60496127381
7 18 2.51027778286
7 19 2.6080036007
8 0 2.58101266264
8 1 2.423807375
8 2 6.85231934507
8 3 8.88026842867
8 4 3.30103580362
8 5 4.77479377435
8 6 7.95781803678
8 7 3.76860034777
8 8 9.54355920945
8 9 3.57568839248
8 10 0.65075560854
8 11 1.01327876912
8 12 4.17245956721
8 13 2.88266379945
8 14 1.57960099252
8 15 8.98917800859
8 16 9.66606961985
8 17 2.12562585429
8 18 4.08416530947
8 19 9.75925690654
9 0 1.07478444585
9 1 3.48715342581
9 2 6.02770114047
9 3 2.56746427249
9 4 5.60273007022
9 5 6.38568228669
9 6 6.61179223593
9 7 9.92671963506
9 8 2.1148493074
9 9 0.910091919462
9 10 0.529179098476
9 11 5.57087640132
9 12 3.67940581918
9 13 6.17797157758
9 14 8.14128471399
9 15 4.87394547038
9 16 3.96815506054
9 17 3.2046833366
9 18 0.560483382065
9 19 5.43547447298
10 0 1.77796986767
10 1 6.94926494541
10 2 9.62285072845
10 3 7.85706785446
10 4 2.53239072339
10 5 6.21037972798
10 6 1.48270023683
10 7 7.98636571324
10 8 9.78104461983
10 9 5.52099891633
10 10 8.41970766737
10 11 7.48138779954
10 12 0.929768946531
10 13 7.48902473765
10 14 1.67683200215
10 15 1.97119417795
10 16 1.13180133036
10 17 4.21291861431
10 18 0.735643931191
10 19 3.92951301782
11 0 9.52495995626
11 1 5.38745842921
11 2 4.12114224081
11 3 9.32067902327
11 4 2.21152599152
11 5 9.95345446887
11 6 7.03052639809
11 7 0.88821662473
11 8 6.93987554638
11 9 4.48330491108
11 10 1.23768384521
11 11 6.07521312221
11 12 0.922499205994
11 13 8.55711207017
11 14 9.71553217245
11 15 6.22672626511
11 16 4.91031287604
11 17 2.51180372654
11 18 3.03542346609
11 19 6.17013816688
12 0 2.5235839251
12 1 7.5593936664
12 2 1.15793105228
12 3 8.87223261638
12 4 8.57323466465
12 5 6.36290188538
12 6 0.34749555651
12 7 9.38171070445
12 8 0.802339065705
12 9 9.03654152746
12 10 7.09911435347
12 11 2.68188603234
12 12 0.481919550751
12 13 5.66292338346
12 14 7.24698445161
12 15 7.81883632797
12 16 7.03192242551
12 17 2.66281630649
12 18 5.6903102191
12 19 5.95610126399
13 0 4.2255605483
13 1 1.5541475274
13 2 2.80448857305
13 3 6.34548317419
13 4 6.77981332859
13 5 4.61550881737
13 6 9.13139478281
13 7 5.62357639458
13 8 4.04270249219
13 9 0.178381828859
13 10 8.15292315012
13 11 6.95995295376
13 12 9.34981083145
13 13 1.50138944331
13 14 6.45066869292
13 15 3.72215699756
13 16 6.0971975762
13 17 2.09039485087
13 18 3.11217136563
13 19 3.45961398968
14 0 2.86360639461
14 1 9.89817741064
14 2 5.72091973204
14 3 0.985399138555
14 4 1.49168114709
14 5 0.605905924771
14 6 7.00393896437
14 7 6.27092620002
14 8 5.24040545073
14 9 3.69500876708
14 10 0.0472974269448
14 11 9.55648196153
14 12 7.19306535913
14 13 3.22069267106
14 14 7.64808390183
14 15 5.36040851866
14 16 3.12308125032
14 17 4.9217742447
14 18 7.41705910308
14 19 7.70104497309
15 0 1.03635697691
15 1 5.16013998747
15 2 6.37987927891
15 3 2.87063563613
15 4 4.38304812645
15 5 6.88547847851
15 6 5.87308529522
15 7 0.691135011962
15 8 7.2450355266
15 9 9.6189105475
15 10 4.1460870474
15 11 0.316176653149
15 12 1.92305182756
15 13 7.01577163902
15 14 6.27203942311
15 15 2.41125265378
15 16 1.74880878107
15 17 2.45286327945
15 18 1.76656468476
15 19 1.05333189499
16 0 9.62992068233
16 1 3.70176648964
16 2 0.188113163549
16 3 8.20007019184
16 4 2.94516162531
16 5 8.42153580027
16 6 7.24924527099
16 7 1.46566937353
16 8 0.10627105338
16 9 8.72927359529
16 10 9.52608957306
16 11 0.755533449082
16 12 2.14505854484
16 13 9.84652106216
16 14 2.66722121468
16 15 1.24493252821
16 16 2.39605197297
16 17 2.94543412501
16 18 9.21314291587
16 19 3.56594338324
17 0 9.93796382799
17 1 3.92056873334
17 2 5.68851583016
17 3 7.24917325949
17 4 7.20654985588
17 5 6.02324070532
17 6 9.8667720524
17 7 8.21201265244
17 8 1.70728478736
17 9 5.48881080074
17 10 0.689156224527
17 11 4.36339802766
17 12 2.48054443734
17 13 7.09681972726
17 14 0.208678813074
17 15 2.9598729733
17 16 1.11552900363
17 17 3.96180053277
17 18 7.67664843153
17 19 8.55499402992
18 0 9.12469375595
18 1 8.71383996386
18 2 3.19036444252
18 3 2.760618474
18 4 5.70984217889
18 5 3.2792607787
18 6 2.3111903649
18 7 0.164159727665
18 8 6.83896115151
18 9 3.27826410775
18 10 7.86770800337
18 11 2.80887262659
18 12 9.52797729781
18 13 3.11720958848
18 14 4.986088299
18 15 3.71775535725
18 16 9.94906937615
18 17 6.08771474961
18 18 2.08024426157
18 19 3.73415403433
19 0 2.85895666174
19 1 5.02896484314
19 2 8.30269764066
19 3 1.09105053967
19 4 3.44614686166
19 5 3.86267452086
19 6 2.83042870284
19 7 1.27706937775
19 8 8.42550234177
19 9 4.86232523132
19 10 3.08247275689
19 11 4.5589056131
19 12 8.35914151704
19 13 8.72728055114
19 14 8.63118303316
19 15 2.08364249169
19 16 4.13329699784
19 17 0.874377228553
19 18 5.02287479277
19 19 1.99505589058
//rootmacro
//non-validated non-validated non-validated non-validated non-validated non-validated
//non-validated non-validated non-validated non-validated non-validated non-validated
//non-validated non-validated non-validated non-validated non-validated non-validated
//original code: https://root.cern.ch/root/html/tutorials/physics/PhaseSpace.C.html
// example of use of TGenPhaseSpace
//Author: Valerio Filippini
//
void DecayEventGen() {
FILE *fp;
char s[256];
TRandom3 *r = new TRandom3();//Mersenne Twistor
double pi = TMath::Pi();
double KE[3],mKE[3];
double theta[3],mtheta[3];
double phi[3],mphi[3];
double dtheta[3] = {1.0, 1.5, 0.5};
double dphi[3] = {0.2, 0.3, 0.1};
double dKE[3] = {0.0002, 0.0003, 0.0004};
TLorentzVector W(0.0, 0.0, 0.0, 3.922565);
//(Momentum, Energy units are Gev/C, GeV)
Double_t masses[3] = { 0.938272, 2.808922, 0.139575};
TGenPhaseSpace event;
event.SetDecay(W, 3, masses);
for (Int_t n=0;n<10;n++) {
Double_t weight = event.Generate();
TLorentzVector *pProton = event.GetDecay(0);
TLorentzVector *pDeuteron = event.GetDecay(1);
TLorentzVector *pPim = event.GetDecay(2);
sprintf(s,"%03d.txt",n);
if ((fp = fopen(s, "w")) == NULL) {
printf("file open error!!\n");
}
for(int p=0; p<3; p++){
TLorentzVector *plv = event.GetDecay(p);
KE[p] = plv->Energy() - masses[p];
mKE[p] = KE[p] + r->Gaus(0.0, dKE[p]);
theta[p] = plv->Theta()/pi*180;
mtheta[p] = theta[p] + r->Gaus(0.0, dtheta[p]);
phi[p] = (pi - plv->Phi())/pi*180;
mphi[p] = phi[p] + r->Gaus(0.0, dphi[p]);
fprintf(fp,"%.4f %.2f\n",mKE[p]*1000, dKE[p]*1000);//in [MeV]
fprintf(fp,"%.4f %.2f\n",mtheta[p], dtheta[p]);
fprintf(fp,"%.4f %.2f\n",mphi[p], dphi[p]);
}
fclose(fp);
}
}
//usage: >root derivative2.c
Double_t myfunc(Double_t *x,Double_t *par)
{
Float_t xx =x[0];
Double_t f = TMath::TanH( par[0]*TMath::Gaus(xx, par[1], par[2]) );
return f;
}
Double_t myder2func(Double_t *x,Double_t *par)
{
TF1* myfuncp = new TF1("myfuncp", myfunc,-10,10,3);
myfuncp->SetParameters(par[0], par[1], par[2]);
return myfuncp->Derivative2(x[0]);
}
void rootmacro_derivative2()
{
double mean = 1.0;
double sigma = 1.5;
double factor =6.0;
TCanvas* c = new TCanvas("c","canvas",800,600);
c->Divide(1,2);
TF1* fp = new TF1("fp",myfunc,-10,10,3);
fp->SetParameters(factor, mean, sigma);
c->cd(1);
fp->Draw();
TF1* fpder2 = new TF1("fpder2",myder2func,-10,10,3);
fpder2->SetParameters(factor, mean, sigma);
c->cd(2);
fpder2->Draw();
c->Print("func.png");
double L3sigma = mean - 3*sigma;
double R3sigma = mean + 3*sigma;
double LMaxPoint = fpder2->GetMaximumX(L3sigma, mean);
double LMinPoint = fpder2->GetMinimumX(L3sigma, mean);
cout << "LMaxPoint: " << LMaxPoint << endl;
cout << "LMinPoint: " << LMinPoint << endl;
double RMinPoint = fpder2->GetMinimumX(mean, R3sigma);
double RMaxPoint = fpder2->GetMaximumX(mean, R3sigma);
cout << "RMinPoint: " << RMinPoint << endl;
cout << "RMaxPoint: " << RMaxPoint << endl;
double LInflectionPoint = fpder2->GetX(0, LMaxPoint, LMinPoint);//second derivative=0
double RInflectionPoint = fpder2->GetX(0, RMinPoint, RMaxPoint);//second derivative=0
double width = RInflectionPoint - LInflectionPoint;
cout << "LInflectionPoint: " << LInflectionPoint << endl;
cout << "RInflectionPoint: " << RInflectionPoint << endl;
cout << "width: " << width << endl;
}
# -*- coding: utf-8 -*-
"""
Created on 20160823
@author: jyoshida-sci
http://y-anz-m.blogspot.jp/2009/07/how-to-use-pyroot-3.html
https://root.cern.ch/phpBB3/viewtopic.php?t=1316
"""
import ROOT
from ROOT import gROOT, TCanvas, TF1
class Linear:
def __call__(self, x, par):
return par[0] + x[0] * par[1]
class myfunc:
def __call__(self, x, par):
return par[0] * ROOT.TMath.TanH( par[1]* ROOT.TMath.Gaus(x[0], par[2], par[3]) )
class myder2func:
def __call__(self, x, par):
myfuncp = TF1("myfuncp", myfunc(),-10,10,4)
myfuncp.SetParameters(par[0], par[1], par[2], par[3]);
return myfuncp.Derivative2(x[0]);
if __name__ == '__main__':
#ROOT
ROOT.std.__file__ = "dummy_for_old_pyastro"
gROOT.Reset()
#Canvas
c1 = TCanvas( 'c1', 'Example with Formula', 200, 10, 700, 500 )
c1.SetGridx()
c1.SetGridy()
#params
factor =1.0
factorg = 0.5
mean = 6.0
sigma = 1.0
#functions
fp = TF1("fp", myfunc(), 0.0, 20.0, 4)
fp.SetParameters(factor, factorg, mean, sigma)
fp.Draw()
fpder2 = TF1("fpder2", myder2func(), 0.0, 20.0, 4)
fpder2.SetParameters(factor, factorg, mean, sigma)
fpder2.Draw("same")
c1.Print('func.png')
#left part
L3sigma = mean - 3*sigma
LMaxPoint = fpder2.GetMaximumX(L3sigma, mean)
LMinPoint = fpder2.GetMinimumX(L3sigma, mean)
print "LMaxPoint: ", LMaxPoint
print "LMinPoint: ", LMinPoint
#right part
R3sigma = mean + 3*sigma
RMinPoint = fpder2.GetMinimumX(mean, R3sigma)
RMaxPoint = fpder2.GetMaximumX(mean, R3sigma)
print "RMinPoint: ", RMinPoint
print "RMaxPoint: ", RMaxPoint
#second derivative=0
LInflectionPoint = fpder2.GetX(0, LMaxPoint, LMinPoint)
RInflectionPoint = fpder2.GetX(0, RMinPoint, RMaxPoint)
print "LInflectionPoint: ", LInflectionPoint
print "RInflectionPoint: ", RInflectionPoint
#width
width = RInflectionPoint - LInflectionPoint
print "width: ", width
//usage: >root this.c
{
const double factor=15;
TArrow* ta = new TArrow(0, 0, 0, 0, 0.005, "|>");
TH2D* h = new TH2D("h","myfield;x;y",50, -1, 11, 50, -1, 12);
gStyle->SetOptStat(0);
h->Draw();
FILE *fp;
char buf[256];
if ((fp = fopen("data.txt", "r")) == NULL) {
printf("file open error!!\n");
return 0;
}
while (fgets(buf, sizeof(buf), fp) != NULL) {
int x,y;
double dx,dy;
sscanf(buf,"%d %d %lf %lf",&x,&y,&dx,&dy);
ta->DrawArrow(x, y, x+dx*factor, y+dy*factor);
}
//reference
ta->DrawArrow(5, 11, 5+0.1*factor, 11);
TText* t=new TText(0, 0,"");
t->SetTextSize(0.035);
t->SetTextColor(kBlack);
t->DrawText(3, 11,"reference");
fclose(fp);
}
//usage: >root this.c
{
TNtuple* nt = new TNtuple("nt","aaaa","a:b:c");
for(int i=0; i<10000; i++){
nt->Fill(1,2,gRandom->Gaus(0,1));
}
TNtuple* nt2 = new TNtuple("nt2","aaaa2","a:b:c");
for(int i=0; i<10000; i++){
nt2->Fill(1,2,gRandom->Gaus(-5.0, 3.0));
}
TH1D *hsum = new TH1D("hsum","histsum",100,-20,20);
nt->Draw("c>>hsum","","");
nt2->Draw("c>>+hsum","","");
}
# -*- coding: utf-8 -*-
"""
Created on Fri Jan 22 18:11:38 2016
"""
import random
for i in range(20):
for j in range(20):
print i, j, random.uniform(0, 10)
//usage: >root this.c
{
for(int i=0; i<100; i++){
x = i*0.1;
y = + 1.0
+ 0.0011*x
+ 0.0022*x*x
+ 0.0033*x*x*x
+ 0.0044*x*x*x*x
+ 0.0055*x*x*x*x*x
+ 0.0066*x*x*x*x*x*x;
printf("%.1f %.5f\n",x,y);
}
}
//usage: >root MeanFreePathOfKminus.c
!!! not validated !!! not validated !!! not validated !!! not validated !!! not validated
/*
Interactions of negative K-mesons in nuclear emulsion at energies (0-180) MeV
B. Bhowmik et al.,
Il Nuovo Cimento Series 10, 16 Febbraio 1964, Volume 31, Issue 4, pp 716-728
DOI:10.1007/BF02733789
Components of the ET-7B emulsion
Quasifree p(K-, K+)Xi- reaction in nuclear emulsion
S. Aoki et al.,
Nuclear Physics A Volume 644, Issue 4, 28 December 1998, Pages 365-385
doi:10.1016/S0375-9474(98)00596-X
*/
//This function returns MeanFreePathOfKminus in the Fuji ET-7B emulsion in cm units
//x: laboratory energy of K-
//par: not used
Double_t MeanFreePathOfKminus_func(Double_t *x,Double_t *par)
{
double Tk = x[0];
double Kmass = 493.7;//[MeV]
double TotalEnergy = Kmass + Tk;
double p = sqrt(TotalEnergy*TotalEnergy - Kmass*Kmass);
double hbar = 197.327 * (10E-13);//[MeV*cm]
double LTk = hbar/p;
//element "H", "C", "N", "O", "S", "Br", "Ag", "I"
double ni[8] = {336.8, 174.0, 49.7, 95.4, 1.4, 93.9, 94.5, 0.5};
double Ai[8] = { 1.0, 12.0, 14.0, 16.0, 32.1, 79.9, 107.9, 126.9};
double Total = 0.0;
for(int i=0 ; i<8; i++){
double r = (1.2E-13)*pow(Ai[i], 1.0/3.0);
double d = r + LTk;
Total += ni[i]*1E20*d*d;
}
return 1.0 / (TMath::Pi()*Total);
}
void MeanFreePathOfKminus()
{
TCanvas* c = new TCanvas("c","canvas",800,600);
double minRange = 0.0;
double maxRange = 2000.0;
int nParams = 0;
TF1* fp = new TF1("fp", MeanFreePathOfKminus_func, minRange, maxRange, nParams);
fp->SetTitle("Mean free path of K^{-};K^{-} laboratory energy [MeV];Mean free path [cm]");
fp->Draw();
double energy = 1000;//[MeV]
printf("f(%.1f[MeV]) = %.2f [cm]\n", energy, MeanFreePathOfKminus_func(&energy,0));
double mfp = 15.0;//[cm]
printf("f(%.1f[MeV]) = %.2f [cm]\n", fp->GetX(mfp), mfp);
}
/*output
f(1600.0[MeV]) = 24.74 [cm]
f(551.9[MeV]) = 15.00 [cm]
*/
{
TH1D* h = new TH1D("h","hist_poisson",20, -0.5, 19.5);
h->SetMinimum(0);
TRandom3* rand = new TRandom3();
for(int i=0; i<1E5; i++){
h->Fill(rand->Poisson(1.5));
}
h->Scale( 1.0 / h->GetEntries());
h->Draw();
TF1 *func = new TF1("func", "TMath::PoissonI(x+0.5,1.5)", 0.0, 10);
func->Draw("same");
TH1D* h2 = new TH1D("h","hist_poisson",10, -0.5, 19.5);
h2->SetMinimum(0);
for(int i=0; i<1E5; i++){
h2->Fill(rand->Poisson(1.5));
}
h2->Scale( 1.0 / h2->GetEntries());
h2->Draw();
vector<Double_t> x,y;
double nu = 1.5;
for(int n=0; n<15; n+=2){
x.push_back( n+0.5 ) ;
y.push_back( TMath::Poisson(n, nu) + TMath::Poisson(n+1, nu) );
}
Double_t* xpointer=&(x.at(0));
Double_t* ypointer=&(y.at(0));
TGraph* tg=new TGraph(x.size(),xpointer,ypointer);
tg->SetMarkerStyle(21);
tg->SetMarkerColor(kRed);
tg->Draw("P same");
//TF1 *func = new TF1("func", "TMath::PoissonI(x+0.5,1.5)", 0.0, 10);
}
0 1
1 4
2 7
3 10
//usage: >root this.c
{
gStyle->SetOptFit(1111);
TGraph* g = new TGraph("pol1.txt","%lg %lg");
g->SetMarkerStyle(8);
g->Draw("AP");//Axis,Point
TF1* f = new TF1("linear_fitting","pol1",0, 25);//"pol1" = First Degree Polynomials
g->Fit("linear_fitting");
f->Draw("same");
}
0.0 1.00000
0.1 1.00014
0.2 1.00034
0.3 1.00067
0.4 1.00120
0.5 1.00206
0.6 1.00347
0.7 1.00574
0.8 1.00931
0.9 1.01482
1.0 1.02310
1.1 1.03526
1.2 1.05271
1.3 1.07724
1.4 1.11109
1.5 1.15696
1.6 1.21815
1.7 1.29859
1.8 1.40295
1.9 1.53670
2.0 1.70620
2.1 1.91883
2.2 2.18304
2.3 2.50848
2.4 2.90614
2.5 3.38837
2.6 3.96912
2.7 4.66396
2.8 5.49026
2.9 6.46733
3.0 7.61650
3.1 8.96134
3.2 10.52775
3.3 12.34412
3.4 14.44152
3.5 16.85382
3.6 19.61787
3.7 22.77369
3.8 26.36460
3.9 30.43746
4.0 35.04280
4.1 40.23505
4.2 46.07269
4.3 52.61851
4.4 59.93977
4.5 68.10839
4.6 77.20122
4.7 87.30024
4.8 98.49274
4.9 110.87160
5.0 124.53550
5.1 139.58916
5.2 156.14357
5.3 174.31627
5.4 194.23156
5.5 216.02078
5.6 239.82258
5.7 265.78314
5.8 294.05651
5.9 324.80482
6.0 358.19860
6.1 394.41705
6.2 433.64832
6.3 476.08983
6.4 521.94855
6.5 571.44130
6.6 624.79507
6.7 682.24734
6.8 744.04637
6.9 810.45155
7.0 881.73370
7.1 958.17545
7.2 1040.07151
7.3 1127.72908
7.4 1221.46814
7.5 1321.62184
7.6 1428.53683
7.7 1542.57364
7.8 1664.10703
7.9 1793.52636
8.0 1931.23600
8.1 2077.65565
8.2 2233.22076
8.3 2398.38295
8.4 2573.61031
8.5 2759.38791
8.6 2956.21813
8.7 3164.62108
8.8 3385.13503
8.9 3618.31683
9.0 3864.74230
9.1 4125.00670
9.2 4399.72514
9.3 4689.53299
9.4 4995.08638
9.5 5317.06261
9.6 5656.16059
9.7 6013.10132
9.8 6388.62835
9.9 6783.50824
//usage: >root this.c
{
gStyle->SetOptFit(1111);
TGraph* g = new TGraph("pol6.txt","%lg %lg");
g->SetMarkerStyle(8);
g->Draw("AP");//Axis,Point
TF1* f = new TF1("linear_fitting","pol6",0, 25);//"pol6" = 6th Degree Polynomials
g->Fit("linear_fitting");
f->Draw("same");
}
#>python pyroot_demo.py
"""
Created on Mon Apr 13 00:28:58 2015
@author: jyoshida
https://root.cern.ch/drupal/content/how-use-use-python-pyroot-interpreter
"""
import ROOT
ROOT.std.__file__ = 'dummy_for_old_pyastro'
if __name__ == '__main__':
ROOT.gROOT.Reset()
c1 = ROOT.TCanvas( 'c1', 'Example with Formula', 200, 10, 700, 500 )
fun1 = ROOT.TF1( 'fun1', 'abs(sin(x)/x)', 0, 10 )
c1.SetGridx()
c1.SetGridy()
fun1.SetNpx(10)#low resolution
fun1.Draw()
c1.Update()
//usage: >root this.c
{
TRandom3 *r = new TRandom3();//Mersenne Twistor
//Gaussian distribution
const double mean = 0.0;
const double sigma = 1.0;
for(int i=0; i<10; i++){
double v = r->Gaus(mean,sigma);
printf("%.4f\n",v);
}
//Uniform distribution
const double maxval = 10.0;
for(int i=0; i<10; i++){
double v = r->Uniform(maxval);
printf("%.4f\n",v);
}
//Uniform Spherical Distribution
const double rad = 10.0;
for(int i=0; i<10; i++){
double x,y,z;
r->Sphere(x, y ,z, rad);
printf("%.4f %.4f %.4f\n",x, y ,z);
}
delete r;
}
//usage: >root this.c
{
vector<double> vd;
vd.push_back(0.0);
vd.push_back(1.0);
vd.push_back(2.0);
vd.push_back(3.0);
vd.push_back(4.0);
//(0+1+2+3+4)/5 = 2.0
double mean = TMath::Mean(vd.size(), &vd.front());
//sqrt(((0-2)^2 + (1-2)^2 + (2-2)^2 + (3-2)^2 + (4-2)^2) / (5-1))
double rms = TMath::RMS(vd.size(), &vd.front());//1.58114
cout << mean << " " << rms << endl;
}
//usage: >root "rootmacro_args.c(""aaaaaa"","1.1","2.2")"
rootmacro_args(const char* str, const double a, const double b){
printf("hello. %s %.1f %.1f\n",str, a, b);
}
//usage: >root this.c
{
TDatabasePDG *pdg = TDatabasePDG::Instance();
TParticlePDG* tp;
printf("name\tMass[GeV]\tLifetime[sec]\tStrangeness\n");
printf("=========================================\n");
tp = pdg->GetParticle(11);
printf("elec\t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(-11);
printf("posi\t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(3122);
printf("Lambda\t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(3312);
printf("Xi- \t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(2212);
printf("p \t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(2112);
printf("n \t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(3112);
printf("Sigma-\t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(111);
printf("pi0 \t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(211);
printf("pi+ \t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
tp = pdg->GetParticle(-211);
printf("pi- \t%e\t%e\t%d\n", tp->Mass(), tp->Lifetime(), tp->Strangeness());
printf("=========================================\n");
printf("Some lifetimes and strangenesses are invalid unless dedicated DB is unloaded.\n");
delete tp;
//do not delete pdg
//delete pdg;
}
//usage: >root this.c
{
TCanvas* c = new TCanvas("c","canvas",800,600);
gStyle->SetOptStat(0);
TH1D *h1 = new TH1D("h1","hist1",50,-10,10);
h1->FillRandom("gaus");//ガウシアン乱数(mean=0.0, sigma=1.0)を5000個
h1->GetXaxis()->SetTitle("#alpha#beta#Gamma#Delta");
h1->GetXaxis()->SetTitleSize(0.07);
h1->Draw();
c->Print("gaus1st.png");
TH1D *h2 = new TH1D("h2","hist2",50,-10,10);
for(int i=0; i<10000; i++){
h2->Fill(gRandom->Gaus(4,0.5));
}
h2->SetTitle("titletitle;xxxxx;yyyyy");
h2->Draw();
c->Print("gaus2nd.png");
TH1D *hsum = new TH1D("hsum","Gaussian distribution;x_title;y_title",50,-10,10);
hsum->Add(h1);
hsum->Add(h2);
hsum->Draw();
c->Print("2gaus.png");
hsum->SetTitle("Normalized Gaussian distribution");
hsum->Scale( 1.0 / hsum->GetEntries());
hsum->Draw();
c->Print("normalized_2gaus.png");
}
//usage: >root this.c
{
//https://root.cern.ch/root/html/TCut.html
TNtuple* nt = new TNtuple("nt","myntuple","aa:bb:cc");
nt->ReadFile("result.txt");
ntr->SetAlias("aaa","aa*1000");
ntr->SetAlias("bbb","bb-20");
TCut mycut = "aaa>10";
ntr->Draw("aaa:cc",mycut,"");
}
//usage: >root this.c
{
TCanvas* c = new TCanvas("c","canvas_title",800,600);
TNtuple* nt = new TNtuple("nt","ntuple_title","a:b");
nt->ReadFile("ab_data.txt");
nt->SetMarkerStyle(6);
nt->Draw("a:b");
c->Print("ab.png");
TH1D* h = new TH1D("h","hist_title",30,-0.5,29.5);
h->SetTitle("Title;XTitle;YTitle");
nt->Draw("a>>h","","");
c->Print("hist.png");
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment