]>
Commit | Line | Data |
---|---|---|
ec5bd527 | 1 | void 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 | ||
18 | void 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 | 117 | void 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 | 271 | void MakePic(){ |
272 | // | |
273 | // | |
274 | // | |
275 | ||
276 | } | |
277 |