streamed selection of the backup track for refit and moved to
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibExB.C
1 void Init(){
2   //
3   // Initialize
4   //
5   AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
6   AliTPCExB::RegisterField(0,field);
7   AliMagF* fieldC0 = new AliMagF("Maps","Maps", 1, 1, AliMagF::k5kG);
8   AliTPCExB::RegisterField(1,fieldC0);
9   AliMagF* fieldC1 = new AliMagF("Maps","Maps", 1, 1, AliMagF::k5kG);
10   AliTPCExB::RegisterField(2,fieldC1);
11
12   gSystem->Load("libSTAT.so");
13   AliTPCExBFirst *exbfirst1  = new  AliTPCExBFirst(fieldC1,0.88*2.6400e+04,50,50,50);
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
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
122   //
123   // Create this file from scan - See AliTPCCalibLaser.C
124   //
125   TFile fscan("laserScan.root");
126   TTree * treeT = (TTree*)fscan.Get("Mean");
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   
139   treeT->SetAlias("side","(-1+(LTr.fP[1]>0)*2)");       // side
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   //
147   fstring+="((dr)^3-1)*bz++";             //1
148   fstring+="((dr)^3-1)*bz*sa++";        //3
149   fstring+="((dr)^3-1)*bz*ca++";        //5
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   //
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   //
159   fstring+="side*((dr)^1-1)*bz++";           //7
160   fstring+="side*((dr)^1-1)*bz*sa++";        //9
161   fstring+="side*((dr)^1-1)*bz*ca++";        //11
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
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());
181   printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
182
183
184   
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)");
187
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)");
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   //
200   TString fstringeb="";
201   //
202   fstringeb+="bz*fbz0++";                //1
203   fstringeb+="sign(bz)*bz^2*fbz2++";     //2
204   fstringeb+="((dr)-1)*bz*sa++";        //9
205   fstringeb+="side*((dr)-1)*bz*sa++";        //9
206
207
208   fstringeb+="((dr)-1)*bz*ca++";        //9
209   fstringeb+="side*((dr)^3-1)*bz*sa++";        //9
210   fstringeb+="side*((dr)^3-1)*bz++";           //7
211
212   fstringeb+="((dr)^3-1)*bz++";           //7
213   fstringeb+="side*((dr)^3-1)*bz*ca++";        //11
214   
215   //  
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);
217   strExB->Tokenize("+")->Print();
218   printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
219   treeT->SetAlias("fitEB",strExB->Data());
220
221
222
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
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  
266 }
267
268
269
270
271 void MakePic(){
272   //
273   //
274   //
275
276 }
277