1 #if !defined(__CINT__) || defined(__MAKECINT__)
5 #include <TLegendEntry.h>
13 #include <TGraphErrors.h>
22 /* Macro to produce LHC file for the luminous region.
23 Author: caffarri@pd.infn.it
27 void VtxInfoLHC(TString vtxtype="TRKnc",
28 TString fname="Vertex.Performance.root",
29 TString ntname="fNtupleBeamSpot"){
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.};
34 Float_t totev=0,totevtriggered=0;
36 Float_t resolVtx=510.;
37 Float_t errResolVtx = 10;
41 Float_t errTotResol =(errResolVtx*errResolVtx)+(0.3*(errp2*errp2)/TMath::Power(meanMult,p2+1));
43 Float_t timeBinning = 120.;
49 Double_t zRange = 40.;
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.;
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);
69 Float_t xVtx,xVtxSc,xerrVtx;
70 Float_t yVtx,yVtxSc,yerrVtx;
71 Float_t zVtx,zVtxSc,zerrVtx;
72 Float_t ntrklets,ncontrVtx,triggered;
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);
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);
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);
98 nt->SetBranchAddress("ntrklets",&ntrklets);
99 trkstitle="Tracklets Multiplicity";
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);
108 Float_t minTstamp=10E+13;
109 Float_t maxTstamp=-1;
112 Float_t minRun = 10E+7;
115 for (Int_t ientries=0; ientries<nt->GetEntriesFast(); ientries++){
116 nt->GetEntry(ientries);
118 if (tstamp<minTstamp) minTstamp=tstamp;
119 if (tstamp>maxTstamp) maxTstamp=tstamp;
121 if (run<minRun) minRun=run-10;
122 if (run>maxRun) maxRun=run+10;
125 Int_t nTimeBins = (Int_t)(maxTstamp-minTstamp) / timeBinning;
127 printf ("number of time bins %d\n", nTimeBins);
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);
134 TH1F **hxVtxMult = new TH1F*[nTimeBins+1]; //
135 TH1F **hyVtxMult = new TH1F*[nTimeBins+1]; //
137 for(Int_t kTimeBin=0;kTimeBin<nTimeBins+1;kTimeBin++){
138 printf ("here we are, kTimeBin = %d\n", kTimeBin);
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);
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);
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);
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);
164 TCanvas *cVtx = new TCanvas("cVtx", "Vtx Canvas", 1000,700);
167 Int_t *fill = new Int_t [nTimeBins+1];
169 for(Int_t iev=0;iev<nnnev;iev++){
172 Int_t rightTimeBin =(Int_t)(tstamp-minTstamp)/ 120;
173 printf ("right Time Bin %d \n", rightTimeBin);
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;
194 // if((tstamp<stopTime)&&(tstamp>=iTstamp)){
196 // ncontrVtx = ntrklets;
198 //printf ("ntrksTRKnc = %d ncontr = %d \n ", ntrklets, ncontrVtx);
200 xVtxSc = 10000.*xVtx;
201 yVtxSc = 10000.*yVtx;
202 zVtxSc = 10000.*zVtx;
205 hxVtx[rightTimeBin]->Fill(xVtx);
206 hyVtx[rightTimeBin]->Fill(yVtx);
207 hzVtx[rightTimeBin]->Fill(zVtx);
209 if((ntrklets>30)&&(ntrklets<45)){
210 hxVtxMult[rightTimeBin]->Fill(xVtx);
211 hyVtxMult[rightTimeBin]->Fill(yVtx);
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;
225 TString outFileName ="aliceDataVertex";
226 TString outFileExt = ".txt";
227 TString outFileFill ;
230 ofstream outFile("1045_lumireg_ALICE.txt");
232 for (Int_t iTimeBin=0; iTimeBin<nTimeBins+1; iTimeBin++){
234 if (fill[iTimeBin] != 1045) continue;
236 Int_t corrTstamp = iTimeBin*120+minTstamp;
238 if (iTimeBin == nTimeBins+1) {
240 hxVtxMult[iTimeBin]->Draw();
242 hyVtxMult[iTimeBin]->Draw();
246 TF1 *fitVtx, *fitVtxMult;
247 if((hxVtx[iTimeBin]->GetEffectiveEntries()<500)){// && (hxVtx[iTimeBin]->GetEffectiveEntries()<500)){
248 //outFile << "0." <<" " << "0."<<" ";
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();
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);
277 if((hyVtx[iTimeBin] ->GetEffectiveEntries()<500)){// && (hyVtxMult[iTimeBin]->GetEffectiveEntries()<500)){
278 //outFile << "0." <<" " << "0."<<" ";
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();
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);
305 if((hzVtx[iTimeBin]->GetEffectiveEntries()<500)){
306 // outFile << "0." <<" " << "0."<<" ";
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();
325 zDiamErr = zSigmaVtxErr;
332 Int_t iBin = htstampX->FindBin(corrTstamp);
333 htstampX->SetBinContent(iBin, xMeanVtx);
334 htstampX->SetBinError(iBin, xMeanVtxErr);
336 iBin = htstampY->FindBin(corrTstamp);
337 htstampY->SetBinContent(iBin, yMeanVtx);
338 htstampY->SetBinError(iBin, yMeanVtxErr);
340 iBin = htstampZ->FindBin(corrTstamp);
341 htstampZ->SetBinContent(iBin, zMeanVtx);
342 htstampZ->SetBinError(iBin, zMeanVtxErr);
344 iBin = htstampSigX->FindBin(corrTstamp);
345 htstampSigX->SetBinContent(iBin, xSigmaVtx);
346 htstampSigX->SetBinError(iBin, xSigmaVtxErr);
348 iBin = htstampSigY->FindBin(corrTstamp);
349 htstampSigY->SetBinContent(iBin, ySigmaVtx);
350 htstampSigY->SetBinError(iBin, ySigmaVtxErr);
352 iBin = htstampSigZ->FindBin(corrTstamp);
353 htstampSigZ->SetBinContent(iBin, zSigmaVtx);
354 htstampSigZ->SetBinError(iBin, zSigmaVtxErr);
356 iBin = htstampXdiam->FindBin(corrTstamp);
357 htstampXdiam->SetBinContent(iBin, xDiam);
358 htstampXdiam->SetBinError(iBin, xDiamErr);
360 iBin = htstampYdiam->FindBin(corrTstamp);
361 htstampYdiam->SetBinContent(iBin,yDiam);
362 htstampYdiam->SetBinError(iBin, yDiamErr);
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;
380 printf ("Err Resol = %f \n", errTotResol);
388 stampCanvasSig->cd(1);
390 stampCanvasSig->cd(2);
392 stampCanvasSig->cd(3);
396 htstampXdiam->Draw();
398 htstampYdiam->Draw();