-
Notifications
You must be signed in to change notification settings - Fork 0
/
myanalysisRemote.C
105 lines (81 loc) · 2.2 KB
/
myanalysisRemote.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
{
TBenchmark bench;
bench.Start("bench");
TChain *t = new TChain("tree");
FILE *f = fopen("listafiles","r");
char nomefile[200];
while(fscanf(f,"%s",nomefile)==1){
t->AddFile(nomefile);
}
// spec: 0=pi+, 1=pi-, 2=K+, 3=K-, 4=p, 5=pbar, ..., 8=phi
// phi -> K+ K-
Int_t ispec1 = 2; // K+
Int_t ispec2 = 3; // K-
Int_t iev=0;
Int_t np=0;
Int_t np1=0;
Int_t np2=0;
Float_t etapart[10000];
Float_t px[10000];
Float_t py[10000];
Float_t pz[10000];
Float_t e[10000];
Int_t id[10000];
Int_t i1list[10000];
Int_t i2list[10000];
Float_t pt,phi,m;
TH2F *h = new TH2F("hPhi","Inv mass distr;p_{T} (GeV/c);m_{KK} (GeV/c^{2});N",20,0,5,1000,0.9,1.9);
for(Int_t i=0;i< t->GetEntries();i++){
t->GetEvent(i);
if(t->GetLeaf("ev")->GetValue() == iev){
// fill particles in the event
etapart[np]=t->GetLeaf("eta")->GetValue();
if(TMath::Abs(etapart[np]) < 1){// && t->GetLeaf("mother")->GetValue() == 8){
pt = t->GetLeaf("pt")->GetValue();
phi = t->GetLeaf("phi")->GetValue();
m = t->GetLeaf("mass")->GetValue();
id[np] = t->GetLeaf("id")->GetValue();
px[np] = pt*TMath::Cos(phi);
py[np] = pt*TMath::Sin(phi);
pz[np] = t->GetLeaf("pz")->GetValue();
e[np] = TMath::Sqrt(m*m + px[np]*px[np] + py[np]*py[np] + pz[np]*pz[np]);
if(id[np] == ispec1){
i1list[np1]=np;
np1++;
}
if(id[np] == ispec2){
i2list[np2]=np;
np2++;
}
np++;
}
}
else{
// process particles filled
for(Int_t ii=0;ii < np1;ii++){
for(Int_t jj=0;jj < np2;jj++){
pt = (px[i1list[ii]] + px[i2list[jj]]) * (px[i1list[ii]] + px[i2list[jj]]);
pt += (py[i1list[ii]] + py[i2list[jj]]) * (py[i1list[ii]] + py[i2list[jj]]);
m = (e[i1list[ii]] + e[i2list[jj]]) * (e[i1list[ii]] + e[i2list[jj]]);
m -= (pz[i1list[ii]] + pz[i2list[jj]]) * (pz[i1list[ii]] + pz[i2list[jj]]);
m -= pt;
pt = TMath::Sqrt(pt);
m = TMath::Sqrt(m);
h->Fill(pt,m);
}
}
// prepare to analyze next event
iev = t->GetLeaf("ev")->GetValue();
i--;
np=0;
np1=0;
np2=0;
}
}
h->Draw("colz");
TFile *fout = new TFile("results.root","RECREATE");
h->Write();
fout->Close();
bench.Stop("bench");
bench.Print("bench");
}