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