Attachment 'lifetime.c'
Download 1 #include<TStyle.h>
2 #include<TFile.h>
3 #include<TTree.h>
4 #include<THStack.h>
5 #include<math.h>
6 #include<TH1.h>
7 #include<RooFit.h>
8 #include<RooPlot.h>
9 #include<RooAbsData.h>
10 #include<RooDataHist.h>
11 #include<RooDataSet.h>
12 #include<RooRealVar.h>
13 #include<RooGaussModel.h>
14 #include<RooDecay.h>
15
16 /******************************************************************/
17 /* lifetime_rooFit calculates the neutralino lifetime, */
18 /* using the reconstructed photonshower and kinematic constraints */
19 /* */
20 /* The fit, a gaussian, convoluted with an exponential */
21 /* is realized with rooFit */
22 /******************************************************************/
23
24 const float c = 299792458; //speed of light [m/s]
25 const float E_cms = 500.0; //center of mass energy [GeV]
26 const float E_0 = 250.0; //beam energy [GeV]
27 const float M_neu1 = 151.0; //neutralino mass [GeV]
28
29 const int MAXFILES = 1; //maximum # cluster per event
30 const int MAXCLUSTER = 30; //maximum # cluster per event
31
32 using namespace RooFit;
33
34 void lifetime_rooFit2() {
35
36 float E_neu1 = E_cms/2; //neutralino energy [GeV]
37 float P_0 = sqrt( pow(E_0,2) - pow(M_neu1,2) ); //neutralino momentum [GeV]
38 float betaGamma = sqrt( pow(E_neu1,2)/pow(M_neu1,2) - 1 ); //beta*gamma
39
40 // give input root tree
41 TFile *myFile[MAXFILES];
42 TString fileName[MAXFILES] = {
43 /* "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_neu1neu1_tau11ns.root"
44 "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n1n1aa.root",
45 "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n2n2aa.root",
46 "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n3n3aa.root",
47 "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_aa.root"*/
48 "/data/nwattime/LDC_optimization/test.root"
49 };
50
51 float nExpected[MAXFILES] = {100};//{2001406,70319,70674,3113697};//{216728};//,2001406,70319,70674,3113697};
52 float efficiency[MAXFILES] ={1};//{1,1,1,1};// {0.5948};//{0.0143,0.0234,0.0233,0.0004};
53 //0.2ns: 0.5948 //2ns: 0.5368 //11ns: 0.1018
54
55 float _weight; float _time;
56
57 // helper tree
58 TFile* _helperFile = new TFile("helperFile","recreate");
59 TTree* _helperTree = new TTree("helperTree","helperTree");
60 _helperTree->SetAutoSave();
61 _helperTree->Branch("time",&_time,"time/F");
62 _helperTree->Branch("weight",&_weight,"weight/F");
63
64 //DEBUG: plot psi
65 TH1F *psiHisto = new TH1F("psi","psi",100,0,600);
66
67 for (int fileCounter=0; fileCounter<MAXFILES; fileCounter++) {
68 myFile[fileCounter] = new TFile(fileName[fileCounter],"READ");
69
70 if ( myFile[fileCounter]->IsOpen() ) {
71
72 int nCluster;
73 float _Egamma[MAXCLUSTER];
74 float _phi[MAXCLUSTER]; float _psi[MAXCLUSTER];
75 float _pos[MAXCLUSTER][3]; float _axis[MAXCLUSTER][3];
76 float _cog[MAXCLUSTER]; float _dir[MAXCLUSTER];
77 float _dist[MAXCLUSTER];
78
79 // open input tree
80 TTree* _tree = (TTree*)myFile[fileCounter]->Get("tree");
81 _tree->SetBranchAddress("cluster_nCluster", &nCluster);
82 _tree->SetBranchAddress("cluster_energy", &_Egamma); //photon energy [GeV]
83 _tree->SetBranchAddress("cluster_clusterAxis", &_axis); //principle axis of inertia [mm]
84 _tree->SetBranchAddress("cluster_position", &_pos); //centre of gravity [mm]
85
86 float nEntries = _tree->GetEntriesFast();
87 _weight = nExpected[fileCounter]*efficiency[fileCounter]/nEntries;
88 std::cout<<"for file: "<<fileName[fileCounter]<<" found weight "<<_weight<<std::endl;
89 for (int entryCounter = 0; entryCounter < nEntries; ++entryCounter) {
90 _tree->GetEntry(entryCounter);
91
92 std::cout<<"Found "<<nCluster<<" clusters in event: "<<entryCounter<<std::endl;
93
94 // calculate reconstructed lifetime
95 for (int clusterCounter=0; clusterCounter<nCluster; ++clusterCounter) {
96
97 // |centre of gravity|
98 _cog[clusterCounter] = sqrt( pow(_pos[clusterCounter][0],2)
99 + pow(_pos[clusterCounter][1],2)
100 + pow(_pos[clusterCounter][2],2) );
101 // |principle axis|
102 _dir[clusterCounter] = sqrt( pow(_axis[clusterCounter][0],2)
103 + pow(_axis[clusterCounter][1],2)
104 + pow(_axis[clusterCounter][2],2) );
105 std::cout<<"cluster cog: "<<_cog[clusterCounter]<<endl;
106 // <(axis,cog) -> cos(psi)= (axis*cog) / (|axis|*|cog|)
107 _psi[clusterCounter] = acos( ( _pos[clusterCounter][0] * _axis[clusterCounter][0]
108 + _pos[clusterCounter][1] * _axis[clusterCounter][1]
109 + _pos[clusterCounter][2] * _axis[clusterCounter][2] )
110 / ( _cog[clusterCounter] * _dir[clusterCounter] ) );
111 // <(neu1,gamma) -> cos(phi) = E0/p0 - m_neu1^2/(2*p0*E_gamma)
112
113 if ( _Egamma[clusterCounter]>0 )
114 _phi[clusterCounter] = acos( E_0/P_0
115 - pow(M_neu1,2) / (2*P_0*_Egamma[clusterCounter]) );
116 else
117 _phi[clusterCounter] = 0;
118
119 //neutralino flight legth lamda = |cog|*sin(psi)/sin(phi)
120 if ( sin(_phi[clusterCounter]) !=0 )
121 _dist[clusterCounter] = _cog[clusterCounter]
122 * sin(_psi[clusterCounter]) / sin(_phi[clusterCounter]); //[mm]
123 else
124 _dist[clusterCounter] = 0;
125
126 _time = _dist[clusterCounter] * 1e6 / (c*betaGamma);
127 _helperTree->Fill();
128
129 std::cout<<"flight dist: "<<_dist[clusterCounter]<<endl;
130 std::cout<<"flight time: "<<_time<<endl;
131
132 //DEBUG: plot psi
133 psiHisto->Fill(_dist[clusterCounter]);
134 }//all reco cluster
135 } //all entries
136 }//file open
137 }//all files
138
139 //DEBUG: plot psi
140 psiHisto->Draw();
141
142 TFile *tree_file = _helperTree->GetCurrentFile();
143 tree_file->Write();
144
145 /********* draw neutralino lifetime spectrum *********/
146
147 RooRealVar time("time","decay time [ns]",0,10);
148 RooRealVar weight("weight","event weight",0,250);
149 // RooDataSet *data = new RooDataSet("data","data",_helperTree,RooArgSet(time,weight),WeightVar(weight));
150 RooDataSet *data = new RooDataSet("data","data",RooArgSet(time,weight),0,WeightVar(weight));
151 // data->setWeightVar(weight);
152 data->Print("v");
153
154 RooPlot* frame = time.frame();
155 data->plotOn(frame,DataError(RooAbsData::SumW2));
156
157 frame->SetTitle("");
158 frame->GetYaxis()->SetTitle("# #chi^{0}_{1}");
159 /*
160 //Decay model (Exponential convoluted with Gauss)
161 //name,title,start value,minimum,maximum
162 RooRealVar mean("mean","resolution mean", 0.270, 0.270, 0.270); //.2,0.,5.
163 RooRealVar sigma("sigma","resolution sigma", 0.185, 0.185, 0.185); //.1,0.,5.
164 RooRealVar decayConst("decayConst","decay constant", 1.0, 0.01, 15.0); //1.,0.,15.
165 RooRealVar fsig("fsig","signal fraction", 0.5, 0.1, 1.0); //0.5,0.,1.
166
167 //for an observable time with uncertainty dt for each event
168 //observable,mean,sigma,observable error
169 RooGaussModel res("res","det. resol.",time,mean,sigma);
170 RooDecay decay("decay","decay model",time,decayConst,res,RooDecay::SingleSided);
171
172 //Gauss model for background: same gaussian as detector resolution
173 //observable,mean,sigma,observable error
174 RooGaussModel BgModel("BgModel","background model.",time,mean,sigma);
175
176 //sum signal and background model
177 //model(x) = fsig*decay(x) + (1-fsig)*BgModel(x)
178 RooAddPdf model("model","model",RooArgList(decay,BgModel),fsig);
179 RooFitResult* r = model->fitTo( *data,
180 Minos(false),
181 Range(0.0, 5.0),
182 Save(true));
183
184 model->plotOn(frame);*/
185 frame->Draw();
186
187 }//lifetime2
188
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.