TPC/CalibMacros/driftITSTPC.C

   1 /*
   2   Formulas:
   3
   4   z  = s* (z0 - vd*(t-t0))
   5
   6   s  - side -1 and +1
   7   t0 - time 0
   8   vd - nominal drift velocity
   9   zs - miscalibrated position
  10
  11   zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
  12   vr  - relative change of the drift velocity
  13   dzt - vd*dt
  14   dr  = zz0-s*z
  15   ..
  16   ==>
  17   zs ~ z - s*vr*(z0-s*z)+s*dzt
  18   --------------------------------
  19   1. Correction function vr constant:
  20
  21
  22   dz = zs-z = -s*vr *(z0-s*z)+s*dzt
  23   dzs/dl = dz/dl +s*s*vr*dz/dl
  24   d(dz/dl) = vr*dz/dl
  25
  26 */
  27
  28
  29 void Init(){
  30
  31   gSystem->Load("libSTAT");
  32   TStatToolkit toolkit;
  33   Double_t chi2;
  34   TVectorD fitParam, fitParam1,fitParam2;
  35   TVectorD errParam, errParam1,errParam2;
  36   TMatrixD covMatrix,covMatrix1,covMatrix2;
  37
  38   Int_t npoints;
  39   //
  40   TFile f("driftitsTPC.root");
  41   TTree * tree = (TTree*)f.Get("Test");
  42   tree->SetAlias("side","(-1+(pTPC.fP[1]>0)*2)");     //side
  43   tree->SetAlias("z","(pTPC.fP[1]+0.0)");     //z position
  44   tree->SetAlias("dr","(1-abs(pTPC.fP[1])/250.)");    //norm drift length
  45   tree->SetAlias("tl","(pTPC.fP[3]+pITS.fP[3])*0.5"); //tan lampbda
  46   tree->SetAlias("sa","sin(pTPC.fAlpha+0.)");         //sin alpha
  47   tree->SetAlias("ca","cos(pTPC.fAlpha+0.)");         //cos alpha
  48
  49   tree->SetAlias("dz","(pTPC.fP[1]-pITS.fP[1])");      //z delta
  50   tree->SetAlias("dy","(pTPC.fP[0]-pITS.fP[0])");      //z delta
  51   tree->SetAlias("dtl","(pTPC.fP[3]-pITS.fP[3])");     //delta tan lampbda
  52   tree->SetAlias("etl","sqrt(pITS.fC[9]+0.)");         //error tan lampbda
  53   tree->SetAlias("ez","sqrt(pITS.fC[2]+0.)");          //error z
  54 }
  55
  56
  57
  58 void InitCuts(){
  59   TCut cut0("abs(dy)<1&&abs(pTPC.fP[1])>10");
  60   //TCut cutOut="(side<0)"
  61   TCut cutA=cut0;
  62   TCut cut38009("run==38009");  // 524
  63   TCut cut38010("run==38010");  // 60
  64   TCut cut38286("run==38286");
  65   TCut cut38532("run==38532");
  66   TCut cut38576("run==38576");
  67   TCut cut38586("run==38586");
  68   TCut cut38591("run==38591");
  69   TCut cut45639("run==45639");
  70   TCut cut46993("run==46993");
  71   TCut cut47164("run==47164");
  72   TCut cut47175("run==47175");
  73   TCut cut47180("run==47180");
  74   TCut cut47274("run==47274");
  75   TCut cutA = cut38591+cut0;
  76   TCut cut1 = "1";
  77   TCut cut2 = "1";
  78 }
  79
  80 //
  81 // variable array
  82 //
  83 npoints=0;
  84 TVectorD vside;
  85 TVectorD vz;
  86 TVectorD vdr;
  87 TVectorD vtl;
  88 TVectorD vsa;
  89 TVectorD vca;
  90 TVectorD vdz;
  91 TVectorD vdtl;
  92 TVectorD vez;
  93 TVectorD vetl;
  94
  95 TVectorD errors;
  96
  97
  98 void FillVar(){
  99   //
 100   //
 101   npoints =tree->Draw("side",cutA+cut1+cut2);
 102   vside.ResizeTo(npoints);  vside.SetElements(tree->GetV1());
 103   //
 104   npoints =tree->Draw("z",cutA+cut1+cut2);
 105   vz.ResizeTo(npoints);  vz.SetElements(tree->GetV1());
 106   //
 107   npoints =tree->Draw("dr",cutA+cut1+cut2);
 108   vdr.ResizeTo(npoints);  vdr.SetElements(tree->GetV1());
 109   //
 110   npoints =tree->Draw("tl",cutA+cut1+cut2);
 111   vtl.ResizeTo(npoints);  vtl.SetElements(tree->GetV1());
 112   //
 113   npoints =tree->Draw("sa",cutA+cut1+cut2);
 114   vsa.ResizeTo(npoints);  vsa.SetElements(tree->GetV1());
 115   //
 116   npoints =tree->Draw("ca",cutA+cut1+cut2);
 117   vca.ResizeTo(npoints);  vca.SetElements(tree->GetV1());
 118   //
 119   npoints =tree->Draw("dz",cutA+cut1+cut2);
 120   vdz.ResizeTo(npoints);  vdz.SetElements(tree->GetV1());
 121   //
 122   npoints =tree->Draw("dtl",cutA+cut1+cut2);
 123   vdtl.ResizeTo(npoints);  vdtl.SetElements(tree->GetV1());
 124   //
 125   npoints =tree->Draw("ez",cutA+cut1+cut2);
 126   vez.ResizeTo(npoints);  vez.SetElements(tree->GetV1());
 127   //
 128   npoints =tree->Draw("etl*0.5",cutA+cut1+cut2);
 129   vetl.ResizeTo(npoints);  vetl.SetElements(tree->GetV1());
 130   //
 131   npoints = tree->Draw("run","run>0");
 132   Int_t *np=new Int_t[npoints]; for (Int_t i=0;i<npoints;i++) np[i]=tree->GetV1()[i];
 133   Int_t *nsort=new Int_t[npoints];
 134   Int_t nruns = toolkit.Freq(npoints,np,nsort,kTRUE);
 135 }
 136
 137
 138 void FitI1(){
 139   //
 140   // Independent fit 1
 141   //
 142   TString fstringTl1="";
 143   fstringTl1+="(side)++";
 144   fstringTl1+="(tl)++";
 145   fstringTl1+="side*(tl)++";
 146   //
 147   //
 148   TString fstringZ1="";
 149   fstringZ1+="(side)++";
 150   fstringZ1+="(-side*(250.0-side*z))++";
 151   fstringZ1+="((250.-side*z))++";
 152   // systematic  shift
 153   fstringZ1+="sa++";
 154   fstringZ1+="ca++";
 155   fstringZ1+="side*sa++";
 156   fstringZ1+="side*ca++";
 157   //
 158   TString *strTl1 = toolkit.FitPlane(tree,"dtl:etl",fstringTl1->Data(), cutA+cut1, chi2,npoints,fitParam,covMatrix);
 159   strTl1->Tokenize("+")->Print();
 160   printf("Tl1: Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 161   tree->SetAlias("fitTl1",strTl1->Data());
 162   tree->Draw("dtl-fitTl1",cutA+cut1);
 163   //
 164   TString *strZ1 = toolkit.FitPlane(tree,"dz:ez",fstringZ1->Data(), cutA+cut1, chi2,npoints,fitParam,covMatrix);
 165   strZ1->Tokenize("+")->Print();
 166   printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 167   tree->SetAlias("fitZ1",strZ1->Data());
 168   tree->Draw("dz-fitZ1",cutA+cut1);
 169   cut1 =  "abs(dz-fitZ1)<2  && abs(dtl-fitTl1)<0.02";
 170 }
 171
 172 void Fit2I(){
 173   //
 174   // Independent fit
 175   //
 176   TString fstringTl2="";
 177   fstringTl2+="(side)++";
 178   fstringTl2+="(tl)++";
 179   fstringTl2+="(side*tl)++";
 180   //
 181   fstringTl2+="(sa)++";
 182   fstringTl2+="(side*sa)++";
 183   fstringTl2+="(ca)++";
 184   fstringTl2+="(side*ca)++";
 185   //
 186   fstringTl2+="(sa*tl)++";
 187   fstringTl2+="(side*sa*tl)++";
 188   fstringTl2+="(ca*tl)++";
 189   fstringTl2+="(side*ca*tl)++";
 190   //
 191   //
 192   //
 193   //
 194   TString fstringZ2="";
 195   fstringZ2+="(side)++";
 196   fstringZ2+="(-side*(250.0-side*z))++";
 197   fstringZ2+="((250.-side*z))++";
 198   // systematic  shift
 199   fstringZ2+="sa++";
 200   fstringZ2+="ca++";
 201   fstringZ2+="side*sa++";
 202   fstringZ2+="side*ca++";
 203   //
 204
 205   TString *strTl2 = toolkit.FitPlane(tree,"dtl:etl",fstringTl2->Data(), cutA+cut1+cut2, chi2,npoints,fitParam,covMatrix);
 206   strTl2->Tokenize("+")->Print();
 207   printf("Tl2: Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 208   tree->SetAlias("fitTl2",strTl2->Data());
 209   tree->Draw("dtl-fitTl2",cutA+cut1);
 210   //
 211   TString *strZ2 = toolkit.FitPlane(tree,"dz:ez",fstringZ2->Data(), cutA+cut1+cut2, chi2,npoints,fitParam,covMatrix);
 212   strZ2->Tokenize("+")->Print();
 213   printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 214   tree->SetAlias("fitZ2",strZ2->Data());
 215   tree->Draw("dz-fitZ2",cutA+cut1);
 216   cut2 =  "abs(dz-fitZ2)<2  && abs(dtl-fitTl2)<0.02";
 217 }
 218
 219
 220 void Fit1(){
 221   //
 222   //
 223   //  dz = zs-z = -s*vr *(z0-s*z)+s*dzt
 224   //  dzs/dl = dz/dl +s*s*vr*dz/dl
 225   //  d(dz/dl) = vr*dz/dl
 226   TLinearFitter fitter1(5, "hyp4");
 227   // parameters
 228   // 0 - P3   offset
 229   // 1 - P1   offset
 230   // 2 - vdr  factor
 231   //
 232   // 3 - side offset P3
 233   // 4 - side offset P1
 234
 235   //  dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s
 236   //  d(dz/dl) = vr*dz/dl               = [0] + [2]*dz/dl        + [3]*s
 237
 238   Double_t xxx[100];
 239   for (Int_t i=0;i<npoints;i++){
 240     for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
 241     // P1
 242     xxx[0] = 1;
 243     xxx[1] = -vside[i]*(250.0-vside[i]*vz[i]);
 244     xxx[3] = vside[i];
 245     fitter1.AddPoint(xxx,vdz[i],vez[i]);
 246     for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
 247     xxx[1] = vtl[i];
 248     xxx[2] = vside[i];
 249     fitter1.AddPoint(xxx,vdtl[i],vetl[i]);
 250   }
 251   fitter1.Eval();
 252   fitter1.GetParameters(fitParam1);
 253   fitter1.GetErrors(errParam1);
 254   //  dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s
 255   //  d(dz/dl) = vr*dz/dl               = [0] + [2]*dz/dl        + [3]*s
 256
 257   TString fstrP1f1 ="";
 258   TString fstrP3f1 ="";
 259   fstrP1f1+=fitParam1[1];fstrP1f1+="-side*(250.0-side*z)*(";fstrP1f1+=fitParam1[2];
 260   fstrP1f1+=")+side*("; fstrP1f1+=fitParam1[4]; fstrP1f1+=")";
 261   //
 262   fstrP3f1+=fitParam1[0];fstrP3f1+="+tl*(";fstrP3f1+=fitParam1[2];
 263   fstrP3f1+=")+side*("; fstrP3f1+=fitParam1[3]; fstrP3f1+=")";
 264   //
 265   tree->SetAlias("fP1f1",fstrP1f1->Data());
 266   tree->SetAlias("fP3f1",fstrP3f1->Data());
 267   //
 268 }
 269
 270
 271
 272 void Fit2(){
 273   //
 274   //
 275   //  dz = zs-z = -s*vr *(z0-s*z)+s*dzt
 276   //  dzs/dl = dz/dl +s*s*vr*dz/dl
 277   //  d(dz/dl) = vr*dz/dl
 278   TLinearFitter fitter2(7, "hyp6");
 279   // parameters
 280   // 0 - P3   offset
 281   // 1 - P1   offset
 282   // 2 - vdr  factor
 283   //
 284   // 3 - side offset P3
 285   // 4 - side offset P1
 286   //
 287   // 5 - vdr difference -A side - C side
 288   // 6 - vdr gradient   -sin(alpha)*tl dependence
 289
 290
 291   //  dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s +[5]* (-s*s*(z0-s*z))
 292   //  dz+= [6]*(-s*(z0-s*z))*sa
 293   //
 294   //  d(dz/dl) = vr*dz/dl               = [0] + [2]*dz/dl        + [3]*s +[5]* s*dz/dl
 295   //   d(dz/dl)+= [6]*dz/dl*sa
 296   Double_t xxx[100];
 297   for (Int_t i=0;i<npoints;i++){
 298     for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
 299     // P1
 300     xxx[0] = 1;
 301     xxx[1] = -vside[i]*(250.0-vside[i]*vz[i]);
 302     xxx[3] = vside[i];
 303     xxx[4] = -vside[i]*vside[i]*(250.0-vside[i]*vz[i]);
 304     xxx[5] = -vside[i]*(250.0-vside[i]*vz[i])*vsa[i];
 305     fitter2.AddPoint(xxx,vdz[i],vez[i]);
 306     //P3
 307     for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
 308     xxx[1] = vtl[i];
 309     xxx[2] = vside[i];
 310     xxx[4] = vside[i]*vtl[i];
 311     xxx[5] = vtl[i]*vsa[i];
 312     fitter2.AddPoint(xxx,vdtl[i],vetl[i]);
 313   }
 314   fitter2.Eval();
 315   fitter2.GetParameters(fitParam2);
 316   fitter2.GetErrors(errParam2);
 317   //  dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s +[5]* (-s*s*(z0-s*z))
 318   //  dz+= [6]*(-s*(z0-s*z))*sa
 319   //
 320   //  d(dz/dl) = vr*dz/dl               = [0] + [2]*dz/dl        + [3]*s +[5]* s*dz/dl
 321   //   d(dz/dl)+= [6]*dz/dl*sa
 322
 323
 324   TString fstrP1f2 ="";
 325   TString fstrP3f2 ="";
 326   fstrP1f2+=fitParam2[1];fstrP1f2+="-side*(250.0-side*z)*(";fstrP1f2+=fitParam2[2];
 327   fstrP1f2+=")+side*("; fstrP1f2+=fitParam2[4];
 328   fstrP1f2+=")-side*side*(250.0-side*z)*("; fstrP1f2+=fitParam2[5];
 329   fstrP1f2+=")-side*(250.0-side*z)*sa*tl*(";fstrP1f2+=fitParam2[6];fstrP1f2+=")";
 330   //
 331   fstrP3f2+=fitParam2[0];fstrP3f2+="+tl*(";fstrP3f2+=fitParam2[2];
 332   fstrP3f2+=")+side*("; fstrP3f2+=fitParam2[3];
 333   fstrP3f2+=")+side*tl*("; fstrP3f2+=fitParam2[5];
 334   fstrP3f2+=")*tl*sa*("; fstrP3f2+=fitParam2[6];fstrP3f2+=")";
 335
 336   //
 337   tree->SetAlias("fP1f2",fstrP1f2->Data());
 338   tree->SetAlias("fP3f2",fstrP3f2->Data());
 339   //
 340 }
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363 TString fstringP3="";
 364 fstringP3+="tl++";
 365 fstringP3+="sa++";
 366 fstringP3+="ca++";
 367 fstringP3+="side*tl++";
 368 fstringP3+="side*sa++";
 369 fstringP3+="side*ca++";
 370
 371 TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA+cut38009, chi2,npoints,fitParam,covMatrix);
 372 strP3->Tokenize("+")->Print();
 373 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 374 tree->SetAlias("fitP3",strP3->Data());
 375 tree->Draw("dtl-fitP3",cut0+cut38009);
 376
 377
 378
 379 TString fstringP3="";
 380 //
 381 fstringP3+="side++";
 382 fstringP3+="tl*(run<40000)++";
 383 fstringP3+="tl*(abs(run-45600)<200)++";
 384 fstringP3+="tl*(run>46000)++";
 385 //
 386 fstringP3+="sa++";
 387 fstringP3+="ca++";
 388 fstringP3+="sa*tl++";
 389 fstringP3+="ca*tl++";
 390 //
 391 fstringP3+="side*tl++";
 392 fstringP3+="side*sa++";
 393 fstringP3+="side*ca++";
 394 fstringP3+="side*sa*tl++";
 395 fstringP3+="side*ca*tl++";
 396
 397
 398
 399 TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA, chi2,npoints,fitParam,covMatrix);
 400 strP3->Tokenize("+")->Print();
 401 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 402 tree->SetAlias("fitP3",strP3->Data());
 403 tree->Draw("dtl-fitP3",cut0);
 404
 405
 406
 407 TString fstringP1="";
 408 //
 409 fstringP1+="side++";
 410 fstringP1+="dr*(run<40000)++";
 411 fstringP1+="dr*(abs(run-45600)<200)++";
 412 fstringP1+="dr*(run>46000)++";
 413 fstringP1+="side*dr*(run<40000)++";
 414 fstringP1+="side*dr*(abs(run-45600)<200)++";
 415 fstringP1+="side*dr*(run>46000)++";
 416 //
 417 fstringP1+="tl*(run<40000)++";
 418 fstringP1+="tl*(abs(run-45600)<200)++";
 419 fstringP1+="tl*(run>46000)++";
 420 fstringP1+="side*tl*(run<40000)++";
 421 fstringP1+="side*tl*(abs(run-45600)<200)++";
 422 fstringP1+="side*tl*(run>46000)++";
 423 //
 424 fstringP1+="sa++";
 425 fstringP1+="ca++";
 426 fstringP1+="sa*tl++";
 427 fstringP1+="ca*tl++";
 428 //
 429 fstringP1+="side*sa++";
 430 fstringP1+="side*ca++";
 431 fstringP1+="side*sa*tl++";
 432 fstringP1+="side*ca*tl++";
 433
 434 TCut cuttl="abs(dtl-fitP3)<0.02"
 435
 436 TString *strP1 = toolkit.FitPlane(tree,"dz",fstringP1->Data(),cuttl+cutA, chi2,npoints,fitParam,covMatrix);
 437 strP1->Tokenize("+")->Print();
 438 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 439 tree->SetAlias("fitP1",strP1->Data());
 440 tree->Draw("dz-fitP1",cuttl+cut0);
 441
 442 */