]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/CalibExB.C
Make OCDB Entry
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibExB.C
CommitLineData
ec5bd527 1void Init(){
2 //
3 // Initialize
4 //
4642ac4b 5 AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
ec5bd527 6 AliTPCExB::RegisterField(0,field);
4642ac4b 7 AliMagF* fieldC0 = new AliMagF("Maps","Maps", 1, 1, AliMagF::k5kG);
58388598 8 AliTPCExB::RegisterField(1,fieldC0);
4642ac4b 9 AliMagF* fieldC1 = new AliMagF("Maps","Maps", 1, 1, AliMagF::k5kG);
58388598 10 AliTPCExB::RegisterField(2,fieldC1);
11
ec5bd527 12 gSystem->Load("libSTAT.so");
58388598 13 AliTPCExBFirst *exbfirst1 = new AliTPCExBFirst(fieldC1,0.88*2.6400e+04,50,50,50);
ec5bd527 14 AliTPCExB::SetInstance(exbfirst1);
15
16}
17
18void FitField(){
19 //
20 //
21 //
22 AliTPCExB * exb = AliTPCExB::Instance();
23 exb->TestExB("field.root");
24 TFile f("field.root");
25 TTree * tree = (TTree*)f.Get("positions");
26 //
27 TStatToolkit toolkit;
28 Double_t chi2;
29 TVectorD fitParam;
30 TMatrixD covMatrix;
31 Int_t npoints;
32 //
33 // SetAliases
34 //
35 tree->SetAlias("sa","sin(phi+0.0)");
36 tree->SetAlias("ca","cos(phi+0.0)");
37 tree->SetAlias("sa2","sin(phi*2+0.0)");
38 tree->SetAlias("ca2","cos(phi*2+0.0)");
39 tree->SetAlias("zn","(x2/250.)");
40 tree->SetAlias("rn","(r/250.)");
41
42 TString fstringSym="";
43 //
44 fstringSym+="zn++";
45 fstringSym+="rn++";
46 fstringSym+="zn*rn++";
47 fstringSym+="zn*zn++";
48 fstringSym+="zn*zn*rn++";
49 fstringSym+="zn*rn*rn++";
50 //
51 fstringSym+="sa++";
52 fstringSym+="ca++";
53 fstringSym+="ca2++";
54 fstringSym+="sa2++";
55 fstringSym+="ca*zn++";
56 fstringSym+="sa*zn++";
57 fstringSym+="ca2*zn++";
58 fstringSym+="sa2*zn++";
59 fstringSym+="ca*zn*zn++";
60 fstringSym+="sa*zn*zn++";
61 fstringSym+="ca*zn*rn++";
62 fstringSym+="sa*zn*rn++";
63
64 // BrBz
65 TString *strBrBz = toolkit.FitPlane(tree,"br/bz",fstringSym->Data(), "abs(r)<250&&abs(r)>90", chi2,npoints,fitParam,covMatrix);
66 strBrBz->Tokenize("+")->Print();
67 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
68 tree->SetAlias("fitBrBz",strBrBz->Data());
69 exb->fMatBrBz = new TVectorD(fitParam);
70 // BrfiBz
71 TString *strBrfiBz = toolkit.FitPlane(tree,"brfi/bz",fstringSym->Data(), "abs(r)<250&&abs(r)>90", chi2,npoints,fitParam,covMatrix);
72 strBrfiBz->Tokenize("+")->Print();
73 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
74 tree->SetAlias("fitBrfiBz",strBrfiBz->Data());
75 exb->fMatBrfiBz = new TVectorD(fitParam);
76
77
78 //
79 // BR integral parameterization
80 //
81 TString *strBRI0 = toolkit.FitPlane(tree,"bri",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2>0&&x2<250", chi2,npoints,fitParam,covMatrix);
82 strBRI0->Tokenize("+")->Print();
83 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
84 tree->SetAlias("fitBRI0",strBRI0->Data());
85 exb->fMatBrBzI0 = new TVectorD(fitParam);
86
87 TString *strBRI1 = toolkit.FitPlane(tree,"bri",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2<-10&&x2>-250", chi2,npoints,fitParam,covMatrix);
88 strBRI1->Tokenize("+")->Print();
89 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
90 tree->SetAlias("fitBRI1",strBRI1->Data());
91 exb->fMatBrBzI1 = new TVectorD(fitParam);
92 //
93
94 TString *strBRFII0 = toolkit.FitPlane(tree,"brfii",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2>0&&x2<250", chi2,npoints,fitParam,covMatrix);
95 strBRFII0->Tokenize("+")->Print();
96 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
97 tree->SetAlias("fitBRFII0",strBRFII0->Data());
98 exb->fMatBrfiBzI0 = new TVectorD(fitParam);
99
100 TString *strBRFII1 = toolkit.FitPlane(tree,"brfii",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2<-10&&x2>-240", chi2,npoints,fitParam,covMatrix);
101 strBRFII1->Tokenize("+")->Print();
102 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
103 tree->SetAlias("fitBRFII1",strBRFII1->Data());
104 exb->fMatBrfiBzI1 = new TVectorD(fitParam);
105
106 TFile * fout = new TFile("bfit.root","recreate");
107 exb->Write("bfit");
108 fout->Close();
109
110}
111
112
113
114
115
116
27f5d48d 117void CalibExB(){
118 //
119 // Macro to test ExB correction versus laser B filed scan data
120 //
121 // Before running it be sure the file laserScan.root exist
27f5d48d 122 //
123 // Create this file from scan - See AliTPCCalibLaser.C
124 //
125 TFile fscan("laserScan.root");
126 TTree * treeT = (TTree*)fscan.Get("Mean");
27f5d48d 127 TStatToolkit toolkit;
128 Double_t chi2;
129 TVectorD fitParam;
130 TMatrixD covMatrix;
131 Int_t npoints;
132
133 TCut cutF0("abs(gphi1-(pphi0+pphi1*bz))<0.05");
134 TCut cutF1("abs(gphiP1-(pphiP0+pphiP1*bz))<0.0005");
135 TCut cutN("entries>2");
136
137 TCut cutA = cutF0+cutF1+cutN;
138
ec5bd527 139 treeT->SetAlias("side","(-1+(LTr.fP[1]>0)*2)"); // side
27f5d48d 140 treeT->SetAlias("dr","(abs(LTr.fP[1]/250.))");
141 treeT->SetAlias("sa","sin(atan2(lx1+0.0,lx0+0.0))");
142 treeT->SetAlias("ca","cos(atan2(lx1+0.0,lx0+0.0))");
143 treeT->SetAlias("ta","tan(asin(LTr.fP[2]+0.0))");
144
145 TString fstring="";
146 //
5fd9157a 147 fstring+="((dr)^3-1)*bz++"; //1
27f5d48d 148 fstring+="((dr)^3-1)*bz*sa++"; //3
27f5d48d 149 fstring+="((dr)^3-1)*bz*ca++"; //5
5fd9157a 150 //
151 fstring+="((dr)^1-1)*bz++"; //1
152 fstring+="((dr)^1-1)*bz*sa++"; //3
153 fstring+="((dr)^1-1)*bz*ca++"; //5
154 //
ec5bd527 155 fstring+="side*((dr)^3-1)*bz++"; //1
156 fstring+="side*((dr)^3-1)*bz*sa++"; //3
157 fstring+="side*((dr)^3-1)*bz*ca++"; //5
158 //
27f5d48d 159 fstring+="side*((dr)^1-1)*bz++"; //7
27f5d48d 160 fstring+="side*((dr)^1-1)*bz*sa++"; //9
27f5d48d 161 fstring+="side*((dr)^1-1)*bz*ca++"; //11
5fd9157a 162 //
163 fstring+="((dr)^3-1)*bz^2*ta++"; //2
164 fstring+="((dr)^3-1)*bz^2*sa*ta++"; //4
165 fstring+="((dr)^3-1)*bz^2*ca*ta++"; //6
166
167 fstring+="((dr)^1-1)*bz^2*ta++"; //2
168 fstring+="((dr)^1-1)*bz^2*sa*ta++"; //4
169 fstring+="((dr)^1-1)*bz^2*ca*ta++"; //6
170 //
171 fstring+="side*((dr)^1-1)*bz^2*ta++"; //8
172 fstring+="side*((dr)^1-1)*bz^2*sa*ta++"; //10
173 fstring+="side*((dr)^1-1)*bz^2*ca*ta++"; //12
27f5d48d 174
175
176
177
178 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix);
179 strq0->Tokenize("+")->Print();
180 treeT->SetAlias("fit",strq0->Data());
ec5bd527 181 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
182
183
27f5d48d 184
58388598 185 treeT->SetAlias("bir1", "-AliTPCExB::GetBrI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
186 treeT->SetAlias("birfi1","-AliTPCExB::GetBrfiI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
ec5bd527 187
58388598 188 treeT->SetAlias("bir0", "AliTPCExB::GetBrI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
189 treeT->SetAlias("birfi0","AliTPCExB::GetBrfiI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
ec5bd527 190
191 treeT->SetAlias("fbz00", "(bir0+birfi0*ta)");
192 treeT->SetAlias("fbz02", "(birfi0+bir0*ta)");
193 treeT->SetAlias("fbz10", "(bir1+birfi1*ta)");
194 treeT->SetAlias("fbz12", "(birfi1+bir1*ta)");
195 //
196 treeT->SetAlias("fbz0", "((fSide==0)*fbz00+(fSide==1)*fbz10)");
197 treeT->SetAlias("fbz2", "((fSide==0)*fbz02+(fSide==1)*fbz12)");
198
199 //
27f5d48d 200 TString fstringeb="";
201 //
ec5bd527 202 fstringeb+="bz*fbz0++"; //1
203 fstringeb+="sign(bz)*bz^2*fbz2++"; //2
ee59d24d 204 fstringeb+="((dr)-1)*bz*sa++"; //9
205 fstringeb+="side*((dr)-1)*bz*sa++"; //9
206
207
208 fstringeb+="((dr)-1)*bz*ca++"; //9
ec5bd527 209 fstringeb+="side*((dr)^3-1)*bz*sa++"; //9
210 fstringeb+="side*((dr)^3-1)*bz++"; //7
ee59d24d 211
212 fstringeb+="((dr)^3-1)*bz++"; //7
213 fstringeb+="side*((dr)^3-1)*bz*ca++"; //11
27f5d48d 214
ec5bd527 215 //
58388598 216 TString *strExB = toolkit.FitPlane(treeT,"gphi1-pphi0",fstringeb->Data(), "abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA, chi2,npoints,fitParam,covMatrix);
27f5d48d 217 strExB->Tokenize("+")->Print();
ec5bd527 218 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
27f5d48d 219 treeT->SetAlias("fitEB",strExB->Data());
ec5bd527 220
221
222
36070491 223 TString fstringeb="";
224 //
225 fstringeb+="AliTPCExB::GetDrphi(254,atan2(lx1,lx0),LTr.fP[1],bz*10)++"; //1
226 fstringeb+="AliTPCExB::GetDr(254,atan2(lx1,lx0),LTr.fP[1],-bz*10)*ta++"; //2
227 //
228 fstringeb+="side*((dr)^1-1)*bz++"; //9
229 fstringeb+="side*((dr)^1-1)*bz*sa++"; //9
230 //
231 fstringeb+="side*((dr)^1-1)*bz*ta++"; //9
232 fstringeb+="side*((dr)^1-1)*bz*sa*ta++"; //9
233 fstringeb+="side*((dr)^1-1)*bz*ca*ta++"; //9
234
235 TString *strExB = toolkit.FitPlane(treeT,"fit",fstringeb->Data(), "abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA, chi2,npoints,fitParam,covMatrix);
236 strExB->Tokenize("+")->Print();
237 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
238 treeT->SetAlias("fitEB",strExB->Data());
239
ec5bd527 240
241
242
243
244
245
246 {
247 Float_t otcor0 =0.7;
248 Float_t omegatau0 =0.35;
249 Float_t rms0 =10;
250 for (Float_t otcor=0.6; otcor<1.1;otcor+=0.1)
251 for (Float_t ometau=0.4; ometau<0.45;ometau+=0.01){
252 char hname[100];
253 sprintf(hname,"hist_ometau_%f(50,-0.10,0.10)",ometau);
254 char expr[100];
255 sprintf(expr,"gphi1-pphi0-bz*fbz0*%f-sign(bz)*(%f*%f*bz)^2*fbz2>>%s",ometau,ometau,otcor,hname);
256 treeT->Draw(expr,"abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA);
257 printf("Ometau=%f\tCor=%f\tRMS=%f\n",ometau,otcor,treeT->GetHistogram()->GetRMS());
258 if (rms0>treeT->GetHistogram()->GetRMS()){
259 otcor0=otcor;
260 omegatau0=ometau;
261 rms0=treeT->GetHistogram()->GetRMS();
262 }
263 }
264 }
265
27f5d48d 266}
ec5bd527 267
268
269
270
58388598 271void MakePic(){
272 //
273 //
274 //
275
276}
277