]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/macros/VtxInfoLHC.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
handle bash args properly using arrays
[u/mrichter/AliRoot.git] / PWGPP / macros / VtxInfoLHC.C
CommitLineData
30b05a14 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <Riostream.h>
3#include <TCanvas.h>
4#include <TLegend.h>
5#include <TLegendEntry.h>
6#include <TFile.h>
7#include <TLine.h>
8#include <TNtuple.h>
9#include <TStyle.h>
10#include <TString.h>
11#include <TSystem.h>
12#include <TH1F.h>
13#include <TGraphErrors.h>
14#include <TF1.h>
15#include <TH2F.h>
16#include <fstream.h>
17#include <TMath.h>
18#include <TProfile.h>
19
20#endif
21
22/* Macro to produce LHC file for the luminous region.
23 Author: caffarri@pd.infn.it
24*/
25
26
27void VtxInfoLHC(TString vtxtype="TRKnc",
28 TString fname="Vertex.Performance.root",
29 TString ntname="fNtupleBeamSpot"){
30
31 // const Int_t nbinsmult=11;
32 // Float_t limmult[nbinsmult+1]={0.,1.5, 2.5, 5.5,10.5,14.5,24.5,30.5,44.5,60.5,74.5, 99999.};
33
34 Float_t totev=0,totevtriggered=0;
35
36 Float_t resolVtx=510.;
37 Float_t errResolVtx = 10;
38 Float_t p2= 1.3;
39 Float_t errp2 = 0.1;
40 Float_t meanMult=38;
41 Float_t errTotResol =(errResolVtx*errResolVtx)+(0.3*(errp2*errp2)/TMath::Power(meanMult,p2+1));
42
43 Float_t timeBinning = 120.;
44
45
46 Char_t hisnam[20];
47 Double_t xRange = 1.;
48 Int_t xBins = 125.;
49 Double_t zRange = 40.;
50 Int_t zBins = 40.;
51
52 /*
53 Int_t binsHRes=250*0.5;
54 Float_t rangeHRes=10000.*0.5;
55 Float_t rangeGrRms = 5000; // micron
56 Float_t rangeGrPull = 15;
57 Float_t rangeHPull=10.;
58 Float_t rangeHResS = 1.;
59 Float_t rangeHResS = 10000*0.5.;
60 */
61
62 TFile *f=new TFile(fname.Data());
63 TList *cOutput = (TList*)f->Get("cOutputVtxESD");
64 TNtuple *nt=(TNtuple*)cOutput->FindObject(ntname.Data());
65 // TNtuple *nt=(TNtuple*)f->FindObject(ntname.Data());
66 Int_t nnnev=nt->GetEntries();
67 printf("Events = %d\n",nnnev);
68
69 Float_t xVtx,xVtxSc,xerrVtx;
70 Float_t yVtx,yVtxSc,yerrVtx;
71 Float_t zVtx,zVtxSc,zerrVtx;
72 Float_t ntrklets,ncontrVtx,triggered;
73 Float_t run, tstamp;
74
75 TString sxx="x"; sxx.Append(vtxtype.Data());
76 nt->SetBranchAddress(sxx.Data(),&xVtx);
77 TString syy="y"; syy.Append(vtxtype.Data());
78 nt->SetBranchAddress(syy.Data(),&yVtx);
79 TString szz="z"; szz.Append(vtxtype.Data());
80 nt->SetBranchAddress(szz.Data(),&zVtx);
81
82 TString xerr="xerr"; xerr.Append(vtxtype.Data());
83 nt->SetBranchAddress(xerr.Data(),&xerrVtx);
84 TString yerr="yerr"; yerr.Append(vtxtype.Data());
85 nt->SetBranchAddress(yerr.Data(),&yerrVtx);
86 TString zerr="zerr"; zerr.Append(vtxtype.Data());
87 nt->SetBranchAddress(zerr.Data(),&zerrVtx);
88
89
90 TCanvas *stampCanvas = new TCanvas ("stampCanvas", "stampCanvas");
91 stampCanvas->Divide(3);
92 TCanvas *stampCanvasSig = new TCanvas ("stampCanvasSig", "stampCanvasSig");
93 stampCanvasSig->Divide(3);
94 TCanvas *diamCanvas = new TCanvas ("diamCanvas", "diamCanvas");
95 diamCanvas->Divide(2);
96
97 TString trkstitle;
98 nt->SetBranchAddress("ntrklets",&ntrklets);
99 trkstitle="Tracklets Multiplicity";
100
101 TString ntrks="ntrks"; ntrks.Append(vtxtype.Data());
102 nt->SetBranchAddress(ntrks.Data(),&ncontrVtx);
103 // nt->SetBranchAddress(ntrks.Data(),&ntrklets);
104 nt->SetBranchAddress("triggered",&triggered);
105 nt->SetBranchAddress("run", &run);
106 nt->SetBranchAddress("tstamp",&tstamp);
107
108 Float_t minTstamp=10E+13;
109 Float_t maxTstamp=-1;
110 Int_t subevents;
111
112 Float_t minRun = 10E+7;
113 Float_t maxRun = 1;
114
115 for (Int_t ientries=0; ientries<nt->GetEntriesFast(); ientries++){
116 nt->GetEntry(ientries);
117
118 if (tstamp<minTstamp) minTstamp=tstamp;
119 if (tstamp>maxTstamp) maxTstamp=tstamp;
120
121 if (run<minRun) minRun=run-10;
122 if (run>maxRun) maxRun=run+10;
123 }
124
125 Int_t nTimeBins = (Int_t)(maxTstamp-minTstamp) / timeBinning;
126
127 printf ("number of time bins %d\n", nTimeBins);
128 subevents=0;
129
130 TH1F **hxVtx = new TH1F*[nTimeBins+1]; //"hxVtx", "hxVtx", 125,-1,1);
131 TH1F **hyVtx = new TH1F*[nTimeBins+1]; //("hyVtx", "hyVtx", 125,-1,1);
132 TH1F **hzVtx = new TH1F*[nTimeBins+1]; //("hzVtx", "hzVtx", 40,-40,40);
133
134 TH1F **hxVtxMult = new TH1F*[nTimeBins+1]; //
135 TH1F **hyVtxMult = new TH1F*[nTimeBins+1]; //
136
137 for(Int_t kTimeBin=0;kTimeBin<nTimeBins+1;kTimeBin++){
138 printf ("here we are, kTimeBin = %d\n", kTimeBin);
139
140 sprintf(hisnam,"hxVtx%d",kTimeBin);
141 hxVtx[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
142 sprintf(hisnam,"hyVtx%d",kTimeBin);
143 hyVtx[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
144 sprintf(hisnam,"hzVtx%d",kTimeBin);
145 hzVtx[kTimeBin]=new TH1F(hisnam,"",zBins,-zRange,zRange);
146 sprintf(hisnam,"hxVtxMult%d",kTimeBin);
147 hxVtxMult[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
148 sprintf(hisnam,"hyVtxMult%d",kTimeBin);
149 hyVtxMult[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
150
151 }
152
153 TH1F *htstampX = new TH1F("tstamp Vx","tstamp Vx", nTimeBins+1, minTstamp, maxTstamp);
154 TH1F *htstampY = new TH1F("tstamp Vy","tstamp Vy", nTimeBins+1, minTstamp, maxTstamp);
155 TH1F *htstampZ = new TH1F("tstamp Vz","tstamp Vz", nTimeBins+1, minTstamp, maxTstamp);
156
157 TH1F *htstampSigX = new TH1F("tstamp Sig x","tstamp Sig x", nTimeBins+1, minTstamp, maxTstamp);
158 TH1F *htstampSigY = new TH1F("tstamp Sig y","tstamp Sig y", nTimeBins+1, minTstamp, maxTstamp);
159 TH1F *htstampSigZ = new TH1F("tstamp Sig z","tstamp Sig z", nTimeBins+1, minTstamp, maxTstamp);
160
161 TH1F *htstampXdiam = new TH1F("tstamp Diam X","tstamp Diam X", nTimeBins+1, minTstamp, maxTstamp);
162 TH1F *htstampYdiam = new TH1F("tstamp Diam Y","tstamp Diam Y", nTimeBins+1, minTstamp, maxTstamp);
163
164 TCanvas *cVtx = new TCanvas("cVtx", "Vtx Canvas", 1000,700);
165 cVtx->Divide(3);
166
167 Int_t *fill = new Int_t [nTimeBins+1];
168
169 for(Int_t iev=0;iev<nnnev;iev++){
170 nt->GetEvent(iev);
171
172 Int_t rightTimeBin =(Int_t)(tstamp-minTstamp)/ 120;
173 printf ("right Time Bin %d \n", rightTimeBin);
174
175 if ((run>=114783) && (run<=114798)) fill[rightTimeBin]=1005;
176 if ((run>=114916) && (run<=114931)) fill[rightTimeBin]=1013;
177 if ((run>=115186) && (run<=115193)) fill[rightTimeBin]=1019;
178 if ((run>=115310) && (run<=115345)) fill[rightTimeBin]=1022;
179 if ((run>=115393) && (run<=115414)) fill[rightTimeBin]=1023;
180 if ((run>=115514) && (run<=115521)) fill[rightTimeBin]=1026;
181 if ((run>=115880) && (run<=115892)) fill[rightTimeBin]=1031;
182 if ((run>=116079) && (run<=116081)) fill[rightTimeBin]=1033;
183 if ((run>=116102) && (run<=116134)) fill[rightTimeBin]=1034;
184 if ((run>=116197) && (run<=116204)) fill[rightTimeBin]=1035;
185 if ((run>=116287) && (run<=116288)) fill[rightTimeBin]=1038;
186 if ((run>=116351) && (run<=116354)) fill[rightTimeBin]=1042;
187 if ((run>=116401) && (run<=116403)) fill[rightTimeBin]=1044;
188 if ((run>=116559) && (run<=116574)) fill[rightTimeBin]=1045;
189 if ((run>=116609) && (run<=116611)) fill[rightTimeBin]=1046;
190 if ((run>=116640) && (run<=116645)) fill[rightTimeBin]=1047;
191 if ((run>=116681) && (run<=116684)) fill[rightTimeBin]=1049;
192
193
194 // if((tstamp<stopTime)&&(tstamp>=iTstamp)){
195 subevents++;
196 // ncontrVtx = ntrklets;
197
198 //printf ("ntrksTRKnc = %d ncontr = %d \n ", ntrklets, ncontrVtx);
199
200 xVtxSc = 10000.*xVtx;
201 yVtxSc = 10000.*yVtx;
202 zVtxSc = 10000.*zVtx;
203
204 if (ncontrVtx>0){
205 hxVtx[rightTimeBin]->Fill(xVtx);
206 hyVtx[rightTimeBin]->Fill(yVtx);
207 hzVtx[rightTimeBin]->Fill(zVtx);
208
209 if((ntrklets>30)&&(ntrklets<45)){
210 hxVtxMult[rightTimeBin]->Fill(xVtx);
211 hyVtxMult[rightTimeBin]->Fill(yVtx);
212 }
213 }
214
215 }
216
217 Double_t xMeanVtx, xMeanVtxErr, xSigmaVtx, xSigmaVtxErr;
218 Double_t yMeanVtx, yMeanVtxErr, ySigmaVtx, ySigmaVtxErr;
219 Double_t zMeanVtx, zMeanVtxErr, zSigmaVtx, zSigmaVtxErr;
220 Double_t xVtxRms, yVtxRms, zVtxRms;
221 Double_t xVtxMult, yVtxMult;
222 Double_t xDiam, xDiamErr, yDiam, yDiamErr, zDiam, zDiamErr;
223
224 /*
225 TString outFileName ="aliceDataVertex";
226 TString outFileExt = ".txt";
227 TString outFileFill ;
228 */
229
230 ofstream outFile("1045_lumireg_ALICE.txt");
231
232 for (Int_t iTimeBin=0; iTimeBin<nTimeBins+1; iTimeBin++){
233
234 if (fill[iTimeBin] != 1045) continue;
235
236 Int_t corrTstamp = iTimeBin*120+minTstamp;
237
238 if (iTimeBin == nTimeBins+1) {
239 cVtx->cd(1);
240 hxVtxMult[iTimeBin]->Draw();
241 cVtx->cd(2);
242 hyVtxMult[iTimeBin]->Draw();
243
244 }
245
246 TF1 *fitVtx, *fitVtxMult;
247 if((hxVtx[iTimeBin]->GetEffectiveEntries()<500)){// && (hxVtx[iTimeBin]->GetEffectiveEntries()<500)){
248 //outFile << "0." <<" " << "0."<<" ";
249 xMeanVtx =0.;
250 xMeanVtxErr = 0.;
251 xSigmaVtx = 0.;
252 xSigmaVtxErr = 0.;
253 xDiam = 0.;
254 xDiamErr =0.;
255 xVtxRms =0.;
256 }
257 else{
258 hxVtx[iTimeBin]->Fit("gaus", "M");
259 fitVtx = hxVtx[iTimeBin]->GetFunction("gaus");
260 xMeanVtx = fitVtx->GetParameter(1);
261 xMeanVtxErr = fitVtx->GetParError(1);
262 xSigmaVtx = fitVtx->GetParameter(2);
263 xSigmaVtxErr = fitVtx->GetParError(2);
264 xVtxRms = hxVtx[iTimeBin]->GetRMS();
265
266 hxVtxMult[iTimeBin]->Fit("gaus", "M");
267 fitVtxMult= hxVtxMult[iTimeBin]->GetFunction("gaus");
268 xVtxMult = fitVtxMult->GetParameter(2)*10000;
269 xDiam=TMath::Sqrt(xVtxMult*xVtxMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
270 if(xDiam >= 500 ) xDiam=0.;
271 xDiamErr = TMath::Sqrt(fitVtxMult->GetParError(2)*10000+ errTotResol);
272 printf ("xDiam = %f errDiam = %f \n", xDiam, xDiamErr);
273
274 }
275
276
277 if((hyVtx[iTimeBin] ->GetEffectiveEntries()<500)){// && (hyVtxMult[iTimeBin]->GetEffectiveEntries()<500)){
278 //outFile << "0." <<" " << "0."<<" ";
279 yMeanVtx=0.;
280 yMeanVtxErr = 0.;
281 ySigmaVtx = 0.;
282 ySigmaVtxErr = 0.;
283 yDiam = 0;
284 yDiamErr =0;
285 yVtxRms=0.;
286 }
287 else{
288 hyVtx[iTimeBin]->Fit("gaus", "M");
289 fitVtx = hyVtx[iTimeBin]->GetFunction("gaus");
290 yMeanVtx = fitVtx->GetParameter(1);
291 yMeanVtxErr = fitVtx->GetParError(1);
292 ySigmaVtx = fitVtx->GetParameter(2);
293 ySigmaVtxErr = fitVtx->GetParError(2);
294 yVtxRms = hyVtx[iTimeBin]->GetRMS();
295
296 hyVtxMult[iTimeBin]->Fit("gaus", "M");
297 fitVtxMult= hyVtxMult[iTimeBin]->GetFunction("gaus");
298 yVtxMult = fitVtxMult->GetParameter(2)*10000;
299 yDiam=TMath::Sqrt(yVtxMult*yVtxMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
300 if(yDiam >= 500 ) yDiam=0.;
301 yDiamErr = TMath::Sqrt(fitVtxMult->GetParError(2)*10000 + errTotResol);
302 printf("yDiam = %f errDiam= %f \n", yDiam, yDiamErr);
303 }
304
305 if((hzVtx[iTimeBin]->GetEffectiveEntries()<500)){
306 // outFile << "0." <<" " << "0."<<" ";
307 zMeanVtx=0.;
308 zMeanVtxErr = 0.;
309 zSigmaVtx = 0.;
310 zSigmaVtxErr = 0.;
311 zDiam = 0.;
312 zDiamErr =0.;
313 zVtxRms =0.;
314 }
315 else{
316 hzVtx[iTimeBin]->Fit("gaus", "M");
317 fitVtx = hzVtx[iTimeBin]->GetFunction("gaus");
318 zMeanVtx = fitVtx->GetParameter(1);
319 zMeanVtxErr = fitVtx->GetParError(1);
320 zSigmaVtx = fitVtx->GetParameter(2);
321 zSigmaVtxErr = fitVtx->GetParError(2);
322 zVtxRms = hzVtx[iTimeBin]->GetRMS();
323
324 zDiam = zSigmaVtx;
325 zDiamErr = zSigmaVtxErr;
326
327
328 }
329
330 //outfile<<
331
332 Int_t iBin = htstampX->FindBin(corrTstamp);
333 htstampX->SetBinContent(iBin, xMeanVtx);
334 htstampX->SetBinError(iBin, xMeanVtxErr);
335
336 iBin = htstampY->FindBin(corrTstamp);
337 htstampY->SetBinContent(iBin, yMeanVtx);
338 htstampY->SetBinError(iBin, yMeanVtxErr);
339
340 iBin = htstampZ->FindBin(corrTstamp);
341 htstampZ->SetBinContent(iBin, zMeanVtx);
342 htstampZ->SetBinError(iBin, zMeanVtxErr);
343
344 iBin = htstampSigX->FindBin(corrTstamp);
345 htstampSigX->SetBinContent(iBin, xSigmaVtx);
346 htstampSigX->SetBinError(iBin, xSigmaVtxErr);
347
348 iBin = htstampSigY->FindBin(corrTstamp);
349 htstampSigY->SetBinContent(iBin, ySigmaVtx);
350 htstampSigY->SetBinError(iBin, ySigmaVtxErr);
351
352 iBin = htstampSigZ->FindBin(corrTstamp);
353 htstampSigZ->SetBinContent(iBin, zSigmaVtx);
354 htstampSigZ->SetBinError(iBin, zSigmaVtxErr);
355
356 iBin = htstampXdiam->FindBin(corrTstamp);
357 htstampXdiam->SetBinContent(iBin, xDiam);
358 htstampXdiam->SetBinError(iBin, xDiamErr);
359
360 iBin = htstampYdiam->FindBin(corrTstamp);
361 htstampYdiam->SetBinContent(iBin,yDiam);
362 htstampYdiam->SetBinError(iBin, yDiamErr);
363
364 outFile <<"2"<<" "<<fill[iTimeBin]<<" "<<corrTstamp<<" ";
365 outFile <<xMeanVtx*10<<" "<<xMeanVtxErr*10<<" ";
366 outFile <<yMeanVtx*10<<" "<<yMeanVtxErr*10<<" ";
367 outFile <<zMeanVtx*10<<" "<<zMeanVtxErr*10<<" ";
368 outFile <<xSigmaVtx*10<<" "<<xSigmaVtxErr*10<<" ";
369 outFile <<ySigmaVtx*10<<" "<<ySigmaVtxErr*10<<" ";
370 outFile <<zSigmaVtx*10<<" "<<zSigmaVtxErr*10<<" ";
371 outFile <<xDiam*0.001<<" "<<xDiamErr*0.001<<" ";
372 outFile <<yDiam*0.001<<" "<<yDiamErr*0.001<<" ";
373 outFile <<zDiam*10<<" "<<zDiamErr*10<<" ";
374 outFile <<xVtxRms*10 <<" "<<yVtxRms*10<<" "<<zVtxRms*10<<endl;
375
376
377
378
379 }
380 printf ("Err Resol = %f \n", errTotResol);
381 stampCanvas->cd(1);
382 htstampX->Draw();
383 stampCanvas->cd(2);
384 htstampY->Draw();
385 stampCanvas->cd(3);
386 htstampZ->Draw();
387
388 stampCanvasSig->cd(1);
389 htstampSigX->Draw();
390 stampCanvasSig->cd(2);
391 htstampSigY->Draw();
392 stampCanvasSig->cd(3);
393 htstampSigZ->Draw();
394
395 diamCanvas->cd(1);
396 htstampXdiam->Draw();
397 diamCanvas->cd(2);
398 htstampYdiam->Draw();
399}