]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/CalibLaserVscan.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
During simulation: fill STU region w/ non null time sums
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibLaserVscan.C
CommitLineData
1a3dcc20 1
2/*
3 //
5813554a 4 // 0. Make a calibration
5 // 1. Make a laser scan list
6 // e.g in TPC workscape
7 find `pwd`/*/laserMean.root >laserScan.txt
8 // 2. Define a reference data
9 rrunA=84469/; for a in `cat laserScan.txt`; do echo `pwd`/$rrunA/laserMean.root; done >laserScanRefA.txt
10 rrunC=84469; for a in `cat laserScan.txt`; do echo `pwd`/$rrunC/laserMean.root; done >laserScanRefC.txt
11 rrun=84469; for a in `cat laserScan.txt`; do echo `pwd`/$rrun/laserMean.root; done >laserScanRef.txt
12 //
13 //
1a3dcc20 14 .x ~/rootlogon.C
15 gSystem->Load("libANALYSIS");
16 gSystem->Load("libTPCcalib");
17 gSystem->Load("libSTAT.so");
18
19 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
20 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
21 .L $ALICE_ROOT/TPC/CalibMacros/CalibLaserVscan.C+
22
23 AliXRDPROOFtoolkit tool;
24 chain = tool.MakeChainRandom("laserScan.txt","Mean",0,10200);
25 chain->Lookup();
26 chainRef = tool.MakeChain("laserScanRef.txt","Mean",0,10200);
27 chain->AddFriend(chainRef,"R")
28 chainRefA = tool.MakeChain("laserScanRefA.txt","Mean",0,10200);
29 chain->AddFriend(chainRefA,"RA")
30 chainRefC = tool.MakeChain("laserScanRefC.txt","Mean",0,10200);
31 chain->AddFriend(chainRefC,"RC")
32 //
33 // MakeMeanBundle();
34 // SaveResult();
35 //
36 ReadRunSetup();
37 ReadResult();
38
39 MakeAnalysisBeam();
40 TFile fbundle("scanDeltaBeam.root");
41 chain->Draw("mphi:GetValueBundle(id,1)","isOK&&GetValueBundle(id,0)>3&&LTr.fSide==0","")
42*/
43
44#include <fstream>
45#include "TFile.h"
46#include "TMatrixD.h"
47#include "TVectorD.h"
48
49#include "TChain.h"
50#include "TCut.h"
51#include "TH1F.h"
52#include "TObjArray.h"
53#include "TProfile.h"
54#include "TCanvas.h"
55#include "TH1F.h"
56#include "TGraph.h"
57#include "TF1.h"
58#include "TTreeStream.h"
59#include "TLegend.h"
60#include "AliTPCLaserTrack.h"
61#include "AliTPCcalibDB.h"
62#include "TStatToolkit.h"
63
64
65TMatrixD matrixP4Corr(7,2); // P4 correction for P0 - y position and P2 -snp(phi)
66TMatrixD matrixP2RMSIROC(7,5); // rms of second derivative IROC
67TMatrixD matrixP2RMSOROC(7,5); // rms of second derivative OROC
68// indexes - beam rod
69// rod 5 means all rods
70//
71TMatrixD matrixMeanBundle(336,6); // mean position/delta per bundle
72TMatrixD matrixRMSBundle(336,6); // rms of distribution
73// indexes:
74// 0 - number of entries
75// 1 - dy
76// 2 - dphi
77// 3 - dp4
78// 4 - dy0
79// 5 - p0:p4 slope
80//
81map<int,TVectorD*> mapRunVoltage; // run to voltage
82// run# -0 ggA 1- ggC 2- coA 3-coC 4-skA 5-skC
83map<int,int> mapRunVgg; // run to gating grid voltage
84map<int,int> mapRunVskirt; //
85map<int,int> mapRunVcover; //
86TArrayI runlist(10000);
87//
88TChain * chain=0;
89TObjArray apic;
90//
91//
92TCut tA="isOK&&RA.isOK";
93TCut tC="isOK&&RC.isOK";
94TCut cA="eY.fElements<0.01&&RA.eY.fElements<0.01&&X.fElements>10&&RA.X.fElements>10";
95TCut cC="eY.fElements<0.01&&RA.eY.fElements<0.01&&X.fElements>10&&RA.X.fElements>10";
96
97
98
99
100Double_t GetValueBundle(Int_t id, Int_t type){
101 //
102 //
103 return matrixMeanBundle(id,type);
104}
105
106Double_t GetP4Corr(Int_t id, Int_t value){
107 AliTPCLaserTrack *ltrack =(AliTPCLaserTrack *) AliTPCLaserTrack::GetTracks()->At(id);
108 Int_t beam = ltrack->GetBeam();
109 return matrixP4Corr(beam,value);
110}
111
112Double_t GetDyBundle(Int_t id){
113 Double_t dy = matrixMeanBundle(id,4);
114 Double_t dx = matrixMeanBundle(id,5);
115 AliTPCLaserTrack *ltrack =(AliTPCLaserTrack *) AliTPCLaserTrack::GetTracks()->At(id);
116 Double_t p2 = ltrack->GetParameter()[2];
117 Double_t delta = dy+dx*p2;
118 return delta;
119}
120
121Double_t GetRMSBundle(Int_t id, Int_t type){
122 //
123 //
124 return matrixRMSBundle(id,type);
125}
126Int_t GetVoltage(Int_t run, Int_t type){
127 //
128 // Get the voltage
129 // run# -0 ggA 1- ggC 2- coA 3-coC 4-skA 5-skC
130 TVectorD *runVoltage = mapRunVoltage[run];
131 if (!runVoltage) return -1;
132 return (*runVoltage)[type];
133}
134
135
136
137void SaveResult(){
138 TFile f("laserData.root","recreate");
139 matrixMeanBundle.Write("matrixMeanBundle");
140 matrixRMSBundle.Write("matrixRMSBundle");
141 matrixP4Corr.Write("matrixP4Corr");
142 matrixP2RMSIROC.Write("matrixP2RMSIROC");
143 matrixP2RMSOROC.Write("matrixP2RMSOROC");
144 f.Close();
145}
146void ReadResult(){
147 AliTPCLaserTrack::LoadTracks();
148 TFile f("laserData.root");
149 matrixMeanBundle.Read("matrixMeanBundle");
150 matrixRMSBundle.Read("matrixRMSBundle");
151 matrixP4Corr.Read("matrixP4Corr");
152 matrixP2RMSIROC.Read("matrixP2RMSIROC");
153 matrixP2RMSOROC.Read("matrixP2RMSOROC");
154 f.Close();
155}
156
157void ReadRunSetup(){
158 ifstream in;
159 in.open("runSetup.txt");
160 TString objfile;
161 string line;
162 TObjArray *arr = 0;
163 Int_t counter=0;
164
165 while(in.good()) {
166 //in >> objfile;
167 getline(in,line);
168 objfile=line;
169 printf("%s\n",objfile.Data());
170 arr = objfile.Tokenize(" ");
171 if (!arr) continue;
172 if (arr->GetEntries()>=7){
173 Int_t run = atoi(arr->At(0)->GetName());
174 Int_t vcover = atoi(arr->At(0)->GetName());
175 Int_t vskirt = atoi(arr->At(1)->GetName());
176 Int_t vgg = atoi(arr->At(3)->GetName());
177 TVectorD * vsetup = new TVectorD(6);
178 for (Int_t iv=0; iv<6; iv++){
179 (*vsetup)[iv] = atoi(arr->At(iv+1)->GetName());
180 }
181 mapRunVoltage[run]=vsetup;
182 mapRunVgg[run]=vgg;
183 mapRunVskirt[run]=vskirt;
184 mapRunVcover[run]=vcover;
185 runlist[counter]=run;
186 counter++;
187 }
188 delete arr;
189 }
190 runlist[counter]=-1;
191}
192
193
194
195
196
197void MakeAnalysisBeam(){
198 //
199 //
200 //
201 TTreeSRedirector *pcstream = new TTreeSRedirector("scanDeltaBeam.root");
202 TH1* phisP0=new TH1F("hhisP0","hhisP0",100,-5.,5.);
203 TH1* phisP0X=new TH1F("hhisP0X","hhisP0X",100,-5.,5.);
204 TH1* phisP2=new TH1F("hhisP2","hhisP2",100,-6,6);
205 TH1* phisP4=new TH1F("hhisP4","hhisP4",100,-0.05,0.05);
206 // second derivative
207 TH1* phisP2IROC=new TH1F("hhisP2IROC","hhisPIROC",100,-2.,2.);
208 TH1* phisP2OROC=new TH1F("hhisP2OROC","hhisPOROC",100,-2.,2.);
209 //
210 //
211 //
212 for (Int_t iside=0;iside<2;iside++){
213 for (Int_t ibeam=0; ibeam<8; ibeam++){
214 for (Int_t irun=0; irun<100; irun++){
215 Int_t run =runlist[irun];
216 if (run==-1) break;
217 if (run==0) continue;
218 Int_t vskirt = GetVoltage(run,1);
219 Int_t vcover = GetVoltage(run,2);
220 Int_t vgg = GetVoltage(run,0);
221 //
222 TCut cut = Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d&&run==%d", iside, ibeam,run);
223 if (ibeam==7){
224 cut = Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&run==%d", iside,run);
225 }
226 Int_t entries = chain->Draw("gp41-GetValueBundle(id,3)>>hhisP4",cut,"goff");
227 if (entries==0) continue;
228 chain->Draw("10*(mphi-GetValueBundle(id,1)-GetP4Corr(id,0)*gp41)>>hhisP0",cut,"goff");
229 chain->Draw("10*(mphi-GetDyBundle(id)-GetP4Corr(id,0)*gp41)>>hhisP0X",cut,"goff");
230 chain->Draw("1000*(mphiP-GetValueBundle(id,2)-GetP4Corr(id,1)*gp41)>>hhisP2",cut,"goff");
231 //
232 printf("%d\t%d\t%d\t%f\t%f\t%d\n", iside, ibeam,run, phisP0X->GetMean(), phisP0X->GetRMS(), Int_t(phisP0X->GetEntries()));
233 //
234 chain->Draw("10*(mphi-GetValueBundle(id,1)-GetP4Corr(id,0)*gp41)>>hhisP0",cut,"goff");
235 //
236 chain->Draw("mPy2vP2In*100^2>>hhisP2IROC",cut,"goff");
237 chain->Draw("mPy2vP2Out*100^2>>hhisP2OROC",cut,"goff");
238 //
239 //
240 Double_t mp0 = phisP0->GetMean();
241 Double_t mp0X = phisP0X->GetMean();
242 Double_t mp2 = phisP2->GetMean();
243 Double_t mp4 = phisP4->GetMean();
244 Double_t mp2I = phisP2IROC->GetMean();
245 Double_t mp2O = phisP2OROC->GetMean();
246 Double_t sp0 = phisP0->GetRMS();
247 Double_t sp0X = phisP0X->GetRMS();
248 Double_t sp2 = phisP2->GetRMS();
249 Double_t sp4 = phisP4->GetRMS();
250 Double_t sp2I = phisP2IROC->GetRMS();
251 Double_t sp2O = phisP2OROC->GetRMS();
252 (*pcstream)<<"vScanBeam"<<
253 "side="<<iside<<
254 "run="<<run<<
255 "ibeam="<<ibeam<<
256 "vgg="<<vgg<<
257 "vskirt="<<vskirt<<
258 "vcover="<<vcover<<
259 "entries="<<entries<<
260 "mp0="<<mp0<<
261 "mp0X="<<mp0X<<
262 "mp2="<<mp2<<
263 "mp4="<<mp4<<
264 "mp2I="<<mp2I<<
265 "mp2O="<<mp2O<<
266 "sp0="<<sp0<<
267 "sp0X="<<sp0X<<
268 "sp2="<<sp2<<
269 "sp4="<<sp4<<
270 "sp2I="<<sp2I<<
271 "sp2O="<<sp2O<<
272 "\n";
273 }
274 }
275 }
276 delete pcstream;
277}
278
279
280
281
282void MakeMeanBundle(){
283 //
284 //
285 //
286 AliTPCLaserTrack::LoadTracks();
287 AliTPCLaserTrack *ltrack;
288 TF1 * fp1 = 0;
289 TCut ccut;
290 TH1F * phisP0 = 0;
291 TH1F * phisP2 = 0;
292 TH1F * phisP4 = 0;
293 //
294 // get p4 corr
295 //
296 for (Int_t ibeam=0; ibeam<7; ibeam++){
297 Int_t entries0=chain->Draw("mphi:gp41",Form("isOK&&LTr.fBeam==%d",ibeam));
298 TGraph gr0(entries0,chain->GetV2(),chain->GetV1());
299 gr0.Draw();
300 gr0.Fit("pol1","Q","Q");
301 fp1 = gr0.GetFunction("pol1");
302 matrixP4Corr(ibeam,0)=fp1->GetParameter(1);
303 Int_t entries2=chain->Draw("mphiP:gp41",Form("isOK&&LTr.fBeam==%d",ibeam));
304 TGraph gr2(entries2,chain->GetV2(),chain->GetV1());
305 gr2.Draw();
306 gr2.Fit("pol1","Q","Q");
307 fp1 = gr2.GetFunction("pol1");
308 matrixP4Corr(ibeam,1)=fp1->GetParameter(1);
309 }
310 //
311 // get RMS of second derivative
312 //
313 phisP0= new TH1F("hisP0","hisP0",500,-1.,1.);
314 for (Int_t irod=0; irod<5;irod++)
315 for (Int_t ibeam=0; ibeam<7; ibeam++){
316 //
317 TCut cut=Form("isOK&&LTr.fBeam==%d&&LTr.fRod==%d",ibeam,irod);
318 if (irod==5) cut=Form("isOK&&LTr.fBeam==%d",ibeam);
319 chain->Draw("mPy2vP2In*100^2>>hisP0",cut);
320 matrixP2RMSIROC(ibeam,irod) = phisP0->GetRMS();
321 chain->Draw("mPy2vP2Out*100^2>>hisP0",cut);
322 matrixP2RMSOROC(ibeam,irod) = phisP0->GetRMS();
323 }
324 delete phisP0;
325
326 //
327 // Get mean bundle position
328 //
329 for (Int_t id=0; id<336; id+=7){
330 ltrack =(AliTPCLaserTrack *) AliTPCLaserTrack::GetTracks()->At(id);
331 Int_t side = ltrack->GetSide();
332 Int_t rod = ltrack->GetRod();
333 Int_t bundle = ltrack->GetBundle();
334 // Int_t beam = ltrack->GetBeam();
335 ccut = Form("isOK&&LTr.fSide==%d&&LTr.fBundle==%d&&LTr.fRod==%d",side,bundle,rod);
336 phisP0= new TH1F("hisP0","hisP0",500,-0.5,0.5);
337 phisP2= new TH1F("hisP2","hisP2",500,-0.02,0.02);
338 phisP4= new TH1F("hisP4","hisP4",500,-0.2,0.2);
339 Int_t entries =chain->Draw("mphi",Form("isOK&&id==%d",id),"goff");
340 chain->Draw("mphi-GetP4Corr(id,0)*gp41>>hisP0",ccut,"");
341 chain->Draw("mphiP-GetP4Corr(id,1)*gp41>>hisP2",ccut,"");
342 chain->Draw("gp41>>hisP4",ccut,"");
343 Int_t entriesG =chain->Draw("mphi-GetP4Corr(id,0)*gp41:tan(asin(LTr.fP[2]))",ccut,"");
344 if (entriesG<3) continue;
345 TGraph gr(entriesG,chain->GetV2(),chain->GetV1());
346 gr.Draw();
347 gr.Fit("pol1","Q","Q");
348 fp1 = gr.GetFunction("pol1");
349 for (Int_t id2=0;id2<6; id2++){
350 Int_t jd = id+id2;
351 matrixMeanBundle(jd,0) = entries;
352 matrixMeanBundle(jd,1) =phisP0->GetMean();
353 matrixRMSBundle(jd,1) =phisP0->GetRMS();
354 matrixMeanBundle(jd,2) =phisP2->GetMean();
355 matrixRMSBundle(jd,2) =phisP2->GetRMS();
356 matrixMeanBundle(jd,3) =phisP4->GetMean();
357 matrixRMSBundle(jd,3) =phisP4->GetRMS();
358 matrixMeanBundle(jd,4) = fp1->GetParameter(0); // delta y
359 matrixMeanBundle(jd,5) = fp1->GetParameter(1); // delta x
360 }
361 //
362 printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\n",id,entries,matrixMeanBundle(id,1),matrixMeanBundle(id,2),matrixMeanBundle(id,3),matrixMeanBundle(id,4),matrixMeanBundle(id,5));
363 printf("%d\t%d\t%f\t%f\t%f\n",id,entries,matrixRMSBundle(id,1),matrixRMSBundle(id,2),matrixRMSBundle(id,3));
364 delete phisP0;
365 delete phisP2;
366 delete phisP4;
367 }
368}
369
370
371
372
373
374
375
376
377void MakeGraphsdY(){
378 //
379 // Make delta Y pictures from voltage scan
380 //
381 TObjArray *aprofY = new TObjArray(14);
382 for (Int_t ib=0;ib<14;ib++){
383 TProfile *profY = new TProfile("py","py",100,0,150);
384 chain->Draw("10*(mphi-GetValueBundle(id,1)):bz>>py",Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d",ib/7,ib%7),"prof");
385 profY->SetName(Form("#DeltaY side %d beam %d",ib/7,ib%7));
386 profY->SetTitle(Form("#DeltaY side %d beam %d",ib/7,ib%7));
387 aprofY->AddAt(profY,ib);
388 profY->SetMaximum(2.5);
389 profY->SetMinimum(-2.5);
390 profY->SetMarkerColor((ib%7)+1);
391 profY->SetMarkerStyle((ib%7)+22);
392 profY->SetMarkerSize(1.2);
393 profY->SetXTitle("U_{gg} (V)");
394 profY->SetYTitle("#Delta_{y} (mm)");
395 profY->Fit("pol1");
396 }
397 TCanvas *cY = new TCanvas("deltaY","deltaY",900,600);
398 cY->Divide(5,3);
399 cY->Draw();
400 {
401 for (Int_t ib=0;ib<14;ib++){
402 cY->cd(ib+1);
403 aprofY->At(ib)->Draw("p");
404 }
405 }
406 cY->SaveAs("pic/deltaYlaserVscna.eps");
407 cY->SaveAs("pic/deltaYlaserVscna.gif");
408 apic.AddLast(cY);
409}
410
411void MakeGraphsdP2(){
412 //
413 // Make delta Y pictures from voltage scan
414 //
415 TObjArray *aprofP2 = new TObjArray(14);
416 for (Int_t ib=0;ib<14;ib++){
417 TProfile *profP2 = new TProfile("pyP","pyP",100,0,150);
418 chain->Draw("(mphiP-GetValueBundle(id,2)):bz>>pyP",Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d",ib/7,ib%7),"prof");
419 profP2->SetName(Form("#Delta_{#phi} side %d beam %d",ib/7,ib%7));
420 profP2->SetTitle(Form("#Delta_{#phi} side %d beam %d",ib/7,ib%7));
421 aprofP2->AddAt(profP2,ib);
422 profP2->SetMaximum(0.005);
423 profP2->SetMinimum(-0.005);
424 profP2->SetMarkerColor((ib%7)+1);
425 profP2->SetMarkerStyle((ib%7)+22);
426 profP2->SetMarkerSize(1.2);
427 profP2->SetXTitle("U_{gg} (V)");
428 profP2->SetYTitle("#Delta_{#phi}");
429 profP2->Fit("pol1");
430 }
431 TCanvas *cP2 = new TCanvas("deltaP2","deltaP2",900,600);
432 cP2->Divide(5,3);
433 cP2->Draw();
434 {
435 for (Int_t ib=0;ib<14;ib++){
436 cP2->cd(ib+1);
437 aprofP2->At(ib)->Draw("p");
438 }
439 }
440 cP2->SaveAs("pic/deltaP2laserVscna.eps");
441 cP2->SaveAs("pic/deltaP2laserVscna.gif");
442 apic.AddLast(cP2);
443}
444
445
446
447
448void MakeGraphsP4(){
449 //
450 // Make delta Y pictures from voltage scan
451 //
452 TObjArray *aprofP4 = new TObjArray(14);
453 for (Int_t ib=0;ib<14;ib++){
454 TProfile *profP4 = new TProfile("pp4","pp4",100,0,150);
455 chain->Draw("(gp41-GetValueBundle(id,3)):bz>>pp4",Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d",ib/7,ib%7),"prof");
456 profP4->SetName(Form("1/p_{t} side %d beam %d",ib/7,ib%7));
457 profP4->SetTitle(Form("1/p_{t} side %d beam %d",ib/7,ib%7));
458 aprofP4->AddAt(profP4,ib);
459 profP4->SetMaximum(0.025);
460 profP4->SetMinimum(-0.025);
461 profP4->SetMarkerColor((ib%7)+1);
462 profP4->SetMarkerStyle((ib%7)+22);
463 profP4->SetMarkerSize(1.2);
464 profP4->SetXTitle("U_{gg} (V)");
465 profP4->SetYTitle("1/p_{t} (GeV/c)");
466 profP4->Fit("pol1");
467 }
468 TCanvas *cY = new TCanvas("P4","P4",900,600);
469 cY->Divide(5,3);
470 cY->Draw();
471 {
472 for (Int_t ib=0;ib<14;ib++){
473 cY->cd(ib+1);
474 aprofP4->At(ib)->Draw("p");
475 }
476 }
477 cY->SaveAs("pic/m1ptlaserVscna.eps");
478 cY->SaveAs("pic/m1ptlaserVscna.gif");
479 apic.AddLast(cY);
480}
481
482
483void MakePlotsP2GG(TCut ucut){
484 //
485 //
486 //
487 TFile fbundle("scanDeltaBeam.root");
488 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
489 TGraph *graph[4];
490 Int_t mstyle[4]={22,23,24,25};
491 Int_t mcolor[4]={2,3,4,6};
492 Int_t entries=0;
493 entries = treeScan->Draw("10*sp2I/4:vgg","ibeam==7&&side==0"+ucut,"");
494 graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
495 entries = treeScan->Draw("10*sp2O/4:vgg","ibeam==7&&side==0"+ucut,"");
496 graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
497 entries = treeScan->Draw("10*sp2I/4:vgg","ibeam==7&&side==1"+ucut,"");
498 graph[2]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
499 entries = treeScan->Draw("10*sp2O/4:vgg","ibeam==7&&side==1"+ucut,"");
500 graph[3]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
501
502 for (Int_t i=0; i<4; i++){
503 graph[i]->GetXaxis()->SetTitle("U_{gg} (V)");
504 graph[i]->GetYaxis()->SetTitle("Sagita (mm)");
505 graph[i]->SetMinimum(0);
506 graph[i]->SetMaximum(2.5);
507 graph[i]->SetMarkerStyle(mstyle[i]);
508 graph[i]->SetMarkerSize(1);
509 graph[i]->SetMarkerColor(mcolor[i]);
510 if (i==0) graph[i]->Draw("ap");
511 graph[i]->Draw("p");
512 }
513 TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan - Sagita on 50 cm");
514 legend->AddEntry(graph[0],"IROC A side");
515 legend->AddEntry(graph[1],"OROC A side");
516 legend->AddEntry(graph[2],"IROC C side");
517 legend->AddEntry(graph[3],"OROC C side");
518 legend->Draw();
519 gPad->SaveAs("pic/scansagita_Vgg.eps");
520 gPad->SaveAs("pic/scansagita_Vgg.gif");
521}
522
523
524
525
526
527void MakePlotsP2Cover(TCut ucut){
528 //
529 //
530 //
531 TFile fbundle("scanDeltaBeam.root");
532 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
533 TGraph *graph[4];
534 Int_t mstyle[4]={22,23,24,25};
535 Int_t mcolor[4]={2,3,4,6};
536 Int_t entries=0;
537 entries = treeScan->Draw("10*sp2I/4:vcover","ibeam==7&&side==0"+ucut,"");
538 graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
539 entries = treeScan->Draw("10*sp2O/4:vcover","ibeam==7&&side==0"+ucut,"");
540 graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
541 entries = treeScan->Draw("10*sp2I/4:vcover","ibeam==7&&side==1"+ucut,"");
542 graph[2]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
543 entries = treeScan->Draw("10*sp2O/4:vcover","ibeam==7&&side==1"+ucut,"");
544 graph[3]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
545
546 for (Int_t i=0; i<4; i++){
547 graph[i]->GetXaxis()->SetTitle("U_{cover} (V)");
548 graph[i]->GetYaxis()->SetTitle("Sagita (mm)");
549 graph[i]->SetMinimum(0.25);
550 graph[i]->SetMaximum(0.7);
551 graph[i]->SetMarkerStyle(mstyle[i]);
552 graph[i]->SetMarkerSize(1);
553 graph[i]->SetMarkerColor(mcolor[i]);
554 if (i==0) graph[i]->Draw("ap");
555 graph[i]->Draw("p");
556 }
557 TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan - Sagita on 50 cm");
558 legend->AddEntry(graph[0],"IROC A side");
559 legend->AddEntry(graph[1],"OROC A side");
560 legend->AddEntry(graph[2],"IROC C side");
561 legend->AddEntry(graph[3],"OROC C side");
562 legend->Draw();
563 gPad->SaveAs("pic/scansagita_Cover.eps");
564 gPad->SaveAs("pic/scansagita_Cover.gif");
565}
566
567void MakePlotsP2Skirt(TCut ucut){
568 //
569 //
570 //
571 TFile fbundle("scanDeltaBeam.root");
572 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
573 TGraph *graph[4];
574 Int_t mstyle[4]={22,23,24,25};
575 Int_t mcolor[4]={2,3,4,6};
576 Int_t entries=0;
577 entries = treeScan->Draw("10*sp2I/4:vskirt","ibeam==7&&side==0"+ucut,"");
578 graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
579 entries = treeScan->Draw("10*sp2O/4:vskirt","ibeam==7&&side==0"+ucut,"");
580 graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
581 entries = treeScan->Draw("10*sp2I/4:vskirt","ibeam==7&&side==1"+ucut,"");
582 graph[2]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
583 entries = treeScan->Draw("10*sp2O/4:vskirt","ibeam==7&&side==1"+ucut,"");
584 graph[3]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
585
586 for (Int_t i=0; i<4; i++){
587 graph[i]->GetXaxis()->SetTitle("U_{skirt} (V)");
588 graph[i]->GetYaxis()->SetTitle("Sagita (mm)");
589 graph[i]->SetMinimum(0.25);
590 graph[i]->SetMaximum(0.7);
591 graph[i]->SetMarkerStyle(mstyle[i]);
592 graph[i]->SetMarkerSize(1);
593 graph[i]->SetMarkerColor(mcolor[i]);
594 if (i==0) graph[i]->Draw("ap");
595 graph[i]->Draw("p");
596 }
597 TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan - Sagita on 50 cm");
598 legend->AddEntry(graph[0],"IROC A side");
599 legend->AddEntry(graph[1],"OROC A side");
600 legend->AddEntry(graph[2],"IROC C side");
601 legend->AddEntry(graph[3],"OROC C side");
602 legend->Draw();
603 gPad->SaveAs("pic/scansagita_Skirt.eps");
604 gPad->SaveAs("pic/scansagita_Skirt.gif");
605}
606
607
608
609
610void MakePlotsdYGG(){
611 //
612 //
613 //
614 TFile fbundle("scanDeltaBeam.root");
615 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
616 TGraph *graph[4];
617 Int_t mstyle[4]={22,23,24,25};
618 Int_t mcolor[4]={2,3,4,6};
619 Int_t entries=0;
620 entries = treeScan->Draw("sp0X:vgg","ibeam==7&&side==0","");
621 graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
622 entries = treeScan->Draw("sp0X:vgg","ibeam==7&&side==1","");
623 graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
624
625 for (Int_t i=0; i<2; i++){
626 graph[i]->GetXaxis()->SetTitle("U_{gg} (V)");
627 graph[i]->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
628 graph[i]->SetMinimum(0);
629 graph[i]->SetMaximum(1);
630 graph[i]->SetMarkerStyle(mstyle[i]);
631 graph[i]->SetMarkerSize(1);
632 graph[i]->SetMarkerColor(mcolor[i]);
633 if (i==0) graph[i]->Draw("ap");
634 graph[i]->Draw("p");
635 }
636 TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan-#sigma_{r#phi}");
637 legend->AddEntry(graph[0],"A side");
638 legend->AddEntry(graph[1],"C side");
639 legend->Draw();
640 gPad->SaveAs("pic/scandy_Vgg.eps");
641 gPad->SaveAs("pic/scandy_Vgg.gif");
642}
643
644
645void MakePlotsdYCover(TCut ucut){
646 //
647 //
648 //
649 TFile fbundle("scanDeltaBeam.root");
650 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
651 TGraph *graph[4];
652 Int_t mstyle[4]={22,23,24,25};
653 Int_t mcolor[4]={2,3,4,6};
654 Int_t entries=0;
655 entries = treeScan->Draw("sp0X:vcover","ibeam==7&&side==0"+ucut,"");
656 graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
657 entries = treeScan->Draw("sp0X:vcover","ibeam==7&&side==1"+ucut,"");
658 graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
659
660 for (Int_t i=0; i<2; i++){
661 graph[i]->GetXaxis()->SetTitle("U_{cover} (V)");
662 graph[i]->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
663 graph[i]->SetMinimum(0);
664 graph[i]->SetMaximum(0.75);
665 graph[i]->SetMarkerStyle(mstyle[i]);
666 graph[i]->SetMarkerSize(1);
667 graph[i]->SetMarkerColor(mcolor[i]);
668 if (i==0) graph[i]->Draw("ap");
669 graph[i]->Draw("p");
670 }
671 TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan-#sigma_{r#phi}");
672 legend->AddEntry(graph[0],"A side");
673 legend->AddEntry(graph[1],"C side");
674 legend->Draw();
675 gPad->SaveAs("pic/scandy_Vcover.eps");
676 gPad->SaveAs("pic/scandy_Vcover.gif");
677}
678
679void MakePlotsdYSkirt(TCut ucut){
680 //
681 //
682 //
683 TFile fbundle("scanDeltaBeam.root");
684 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
685 TGraph *graph[4];
686 Int_t mstyle[4]={22,23,24,25};
687 Int_t mcolor[4]={2,3,4,6};
688 Int_t entries=0;
689 entries = treeScan->Draw("sp0X:vskirt","ibeam==7&&side==0"+ucut,"");
690 graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
691 entries = treeScan->Draw("sp0X:vskirt","ibeam==7&&side==1"+ucut,"");
692 graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
693
694 for (Int_t i=0; i<2; i++){
695 graph[i]->GetXaxis()->SetTitle("U_{skirt} (V)");
696 graph[i]->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
697 graph[i]->SetMinimum(0);
698 graph[i]->SetMaximum(0.75);
699 graph[i]->SetMarkerStyle(mstyle[i]);
700 graph[i]->SetMarkerSize(1);
701 graph[i]->SetMarkerColor(mcolor[i]);
702 if (i==0) graph[i]->Draw("ap");
703 graph[i]->Draw("p");
704 }
705 TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan-#sigma_{r#phi}");
706 legend->AddEntry(graph[0],"A side");
707 legend->AddEntry(graph[1],"C side");
708 legend->Draw();
709 gPad->SaveAs("pic/scandy_Vskirt.eps");
710 gPad->SaveAs("pic/scandy_Vskirt.gif");
711}
712
713
714
715void GetOptimalSetting(){
716 //
717 //
718 //
719 TFile fbundle("scanDeltaBeam.root");
720 TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
721 //
722 Float_t meansp0X=0;
723 Float_t meansp2I=0;
724 Float_t meansp2O=0;
725 TH1F * his = new TH1F("hisDelta","hisDelta",300,0,1);
726 treeScan->Draw("sp0X>>hisDelta","ibeam==7","");
727 meansp0X = his->GetMean();
728 treeScan->Draw("sp2I>>hisDelta","ibeam==7","");
729 meansp2I = his->GetMean();
730 treeScan->Draw("sp2O>>hisDelta","ibeam==7","");
731 meansp2O = his->GetMean();
732 //
733 //
734 //
735 Int_t indexes[10000];
736 for (Int_t iside=0;iside<2;iside++){
737 printf("*************************************\n");
738 printf("Side%d\tRun\tchi2\t\tVgg\tVskirt\tVcover\n",iside);
739
740 Int_t entries = treeScan->Draw(Form("sp0X/%f+sp2I/%f+sp2O/%f:run",meansp0X,meansp2I,meansp2O),Form("ibeam==7&&side==%d",iside));
741 TMath::Sort(entries,treeScan->GetV1(),indexes,kFALSE);
742 for (Int_t i=0; i<entries; i++){
743 Int_t index = indexes[i];
744 Int_t run = treeScan->GetV2()[index];
745 Float_t chi2=treeScan->GetV1()[index];
746 Float_t vgg =
747 printf("%d\t%d\t%f\t%d\t%d\t%d\n",iside,run, chi2,GetVoltage(run,0), GetVoltage(run,1),GetVoltage(run,2));
748 }
749 }
750
751 for (Int_t iside=0;iside<2;iside++){
752 printf("***********************\n");
753 printf("Side%d\tRun\tchi2\t\tVgg\tVskirt\tVcover\n",iside);
754
755 Int_t entries = treeScan->Draw(Form("sp2I/%f+sp2O/%f:run",meansp0X,meansp2I,meansp2O),Form("ibeam==7&&side==%d",iside));
756 TMath::Sort(entries,treeScan->GetV1(),indexes,kFALSE);
757 for (Int_t i=0; i<entries; i++){
758 Int_t index = indexes[i];
759 Int_t run = treeScan->GetV2()[index];
760 Float_t chi2=treeScan->GetV1()[index];
761 Float_t vgg =
762 printf("%d\t%d\t%f\t%d\t%d\t%d\n",iside,run, chi2,GetVoltage(run,0), GetVoltage(run,1),GetVoltage(run,2));
763 }
764 }
765}
766
767
768
769
770void MakeAliases(){
771 //
772 // use table
773 //
774 chain->SetAlias("VggA","GetVoltage(run,0)");
775 chain->SetAlias("VggC","GetVoltage(run,1)");
776 chain->SetAlias("VcoA","GetVoltage(run,2)");
777 chain->SetAlias("VcoC","GetVoltage(run,3)");
778 chain->SetAlias("VskA","GetVoltage(run,4)");
779 chain->SetAlias("VskC","GetVoltage(run,5)");
780 //
781 // cuts
782 //
783 chain->SetAlias("TisOK","mdEdx>5&&entries>400");
784 chain->SetAlias("CisOK","nCl.fElements>entries*0.5&&eY.fElements<0.01");
785 chain->SetAlias("ATisOK","RA.mdEdx>5&&RA.entries>400");
786 chain->SetAlias("ACisOK","RA.nCl.fElements>RA.entries*0.5&&abs(RA.dY.fElements-dY.fElements)<0.3&&RA.eY.fElements<0.01");
787 //
788 // shortcuts
789 //
790 chain->SetAlias("lX","X.fElements"); //
791 chain->SetAlias("tY","kY.fElements"); //
792 chain->SetAlias("dE","sign(Y.fElements)*(tan(10*pi/180.)*X.fElements-abs(Y.fElements))");
793 chain->SetAlias("RdY","(dY.fElements-RA.dY.fElements)");
794
795 chain->SetAlias("dIF","lX-85.2");
796 chain->SetAlias("dOF","245.8-lX");
797 chain->SetAlias("dUtY","((GetVoltage(run,0)-GetVoltage(RA.run,0))/400.*tY)");
798 //
799 chain->SetAlias("cIF","dUtY/(1+dIF/13.0)");
800 chain->SetAlias("cOF","dUtY/(1+dOF/22.0)");
801 chain->SetAlias("cIF2","(dUtY*dIF/13.0)/((1.+dIF/13.0)^2)");
802 chain->SetAlias("cOF2","(dUtY*dOF/22.0)/((1.+dOF/22.0)^2)");
803
804 chain->SetAlias("cE","(GetVoltage(run,0)-GetVoltage(RA.run,0))/400.*sign(dE)/(1+(abs(dE)-1.55)/1.4)");
805 chain->SetAlias("cE2","(GetVoltage(run,0)-GetVoltage(RA.run,0))/400.*sign(dE)*((abs(dE)-1.55)/1.4)/((1+(abs(dE)-1.55)/1.4)^2)");
806}
807
808
809
810void MakeAliasesBoth(){
811 //
812 // cuts - slect good tracks
813 //
814 chain->SetAlias("TisOK","mdEdx>5&&entries>400");
815 chain->SetAlias("ATisOK","(LTr.fSide==0)*(RA.mdEdx>5&&RA.entries>500)");
816 chain->SetAlias("CTisOK","(LTr.fSide==1)*(RC.mdEdx>5&&RC.entries>500)");
817 //
818 chain->Draw(">>run","TisOK&&(ATisOK||CTisOK)","entryList");
819 TEntryList *elist = (TEntryList*)gDirectory->Get("run");
820 chain->SetEntryList(elist);
821 //
822 //
823 chain->SetAlias("CisOK","(nCl.fElements>entries*0.5&&eY.fElements<0.01)");
824 chain->SetAlias("ACisOK","RA.nCl.fElements>RA.entries*0.5&&abs(RA.dY.fElements-dY.fElements)<0.3&&RA.eY.fElements<0.01");
825 chain->SetAlias("CCisOK","RC.nCl.fElements>RC.entries*0.5&&abs(RC.dY.fElements-dY.fElements)<0.3&&abs(RC.dZ.fElements-dZ.fElements)<0.3&&RC.eY.fElements<0.01");
826
827
828 //
829 // voltage table
830 //
831 chain->SetAlias("Vgg","((LTr.fSide==0)*GetVoltage(run,0)+(LTr.fSide==1)*GetVoltage(run,1))");
832 chain->SetAlias("VcoA","((LTr.fSide==0)*GetVoltage(run,2)+(LTr.fSide==1)*GetVoltage(run,3))");
833 chain->SetAlias("VskA","((LTr.fSide==0)*GetVoltage(run,4)+(LTr.fSide==1)*GetVoltage(run,5))");
834
835 chain->SetAlias("dVgg","(((LTr.fSide==0)*(GetVoltage(run,0)-GetVoltage(RA.run,0))+((LTr.fSide==0)*(GetVoltage(run,0)-GetVoltage(RA.run,0))))/400.)");
836 //
837 // shortcuts
838 //
839 chain->SetAlias("RdY","((LTr.fSide==0)*(dY.fElements-RA.dY.fElements)+(LTr.fSide==1)*(dY.fElements-RC.dY.fElements))");
840
841 chain->SetAlias("drift","(1.-abs(LTr.fP[1]+0.)/250)");
842 chain->SetAlias("lX","X.fElements"); //
843 chain->SetAlias("tY","kY.fElements"); //
844 chain->SetAlias("dE","sign(Y.fElements)*(tan(10*pi/180.)*X.fElements-abs(Y.fElements))");
845
846 chain->SetAlias("dIF","lX-85.2");
847 chain->SetAlias("dOF","245.8-lX");
848 chain->SetAlias("dUtY","dVgg*tY");
849 //
850 //
851 //
852 chain->SetAlias("cIF","dUtY/(1+dIF/13.0)");
853 chain->SetAlias("cOF","dUtY/(1+dOF/22.0)");
854 chain->SetAlias("cIF2","(dUtY*dIF/13.0)/((1.+dIF/13.0)^2)");
855 chain->SetAlias("cOF2","(dUtY*dOF/22.0)/((1.+dOF/22.0)^2)");
856
857 chain->SetAlias("cE","dVgg*sign(dE)/(1+(abs(dE)-1.5)/1.3)");
858 chain->SetAlias("cE2","dVgg*sign(dE)*((abs(dE)-1.5)/1.3)/((1+(abs(dE)-1.5)/1.3)^2)");
859
860}
861
862
863
864
865void MakeFit(){
866
867 Int_t ntracks=3000000;
868 TStatToolkit toolkit;
869 Double_t chi2=0;
870 Int_t npoints=0;
871 TVectorD fitParam;
872 TMatrixD covMatrix;
873 TString fstring="";
874 fstring+="cIF++";
875 fstring+="cIF2++";
876 fstring+="cOF++";
877 fstring+="cOF2++";
878 fstring+="cE++";
879 fstring+="cE2++";
880 TCut cutA="LTr.fBeam>-1&&CisOK&&(ACisOK||CCisOK)";
881
882 TString * strA = TStatToolkit::FitPlane(chain,"RdY", fstring.Data(),cutA, chi2,npoints,fitParam,covMatrix,-1.,0, ntracks);
883 chain->SetAlias("fdY",strA->Data());
884
885 printf("sqrt(Chi2/npoints)=%f\n",TMath::Sqrt(chi2/npoints));
886 chain->Draw("RdY-fdY",cutA);
887
888
889 TF1 fe("fe","[0]/(1+(x-[2])/[1])",2,10);
890 fe.SetParameters(0.4,1,1.5);
891
892
893}
894
895
896
897//Examples:
898// chain->Draw("RdY:(GetVoltage(run,0)-GetVoltage(RA.run,0))*tY/(1+dIF/10.)","TisOK&&ATisOK&&CisOK&&ACisOK&&LTr.fBundle>0&&LTr.fBeam!=3&&dIF<30","",10000)
899
900
901
902