TPC/CalibMacros/AliTPCCorrectionDemo.C

   1 /*
   2  gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TPC -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TOF -I$ALICE_ROOT/RAW -I$ALICE_ROOT/STAT");
   3
   4  */
   5 #if !defined(__CINT__) || defined(__MAKECINT__)
   6 #include "THnSparse.h"
   7 #include "TLatex.h"
   8 #include "TCanvas.h"
   9 #include "TLegend.h"
  10 #include "TSystem.h"
  11 #include "TFile.h"
  12 #include "TChain.h"
  13 #include "TCut.h"
  14 #include "TH3.h"
  15 #include "TH2F.h"
  16 #include "TProfile3D.h"
  17 #include "TMath.h"
  18 #include "TVectorD.h"
  19 #include "TMatrixD.h"
  20 #include "TStatToolkit.h"
  21 #include "TTreeStream.h"
  22 #include "AliExternalTrackParam.h"
  23 #include "AliESDfriend.h"
  24 #include "AliTPCcalibTime.h"
  25 #include "TROOT.h"
  26 #include "AliXRDPROOFtoolkit.h"
  27 #include "AliTPCCorrection.h"
  28 #include "AliTPCExBTwist.h"
  29 #include "AliTPCGGVoltError.h"
  30 #include "AliTPCComposedCorrection.h"
  31 #include "AliTPCExBConical.h"
  32 #include "TPostScript.h"
  33 #include "TStyle.h"
  34 #include "AliTrackerBase.h"
  35 #include "TGraph.h"
  36 #include "AliCDBManager.h"
  37 #include "AliTPCExBBShape.h"
  38 #include "TRandom.h"
  39 #include "AliGeomManager.h"
  40 #include "AliESDVertex.h"
  41 #include "AliTPCcalibDB.h"
  42 #endif
  43
  44
  45 void AliTPCCorrectionDemo() {
  46
  47   //
  48   // This is a Demo function of the general class AliTPCCorrection, which is used for
  49   // general space point correction due to different effects.
  50   // The effects used in this Demo are:
  51   //   1. ExB twist - general offset of the TPC axis in comparison to the B field axis
  52   //   2. GG error (Gating Grid volt. error) - not perfectly aligned GG voltage (in terms of voltage)
  53   //   3. ExBBShape - B field shape correction of the secound order
  54   //
  55   // See class descriptions for further details
  56   //
  57   // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
  58   //
  59   //
  60   // omegaTau (wt) of the langevin equation
  61   // This is a function of the drift vel., the magnetic and electric field
  62   // e.g. vd=2.6 cm/usc; Ez=400 V/cm; Bz=0.5 T
  63   // wt =  -10.0*(Bz*10)*vd/Ez = -0.325
  64
  65   Double_t vdrift = 2.6; // [cm/us]   // to be updated: per second (ideally)
  66   Double_t bzField = -0.5; // [Tesla] // to be updated: per run
  67   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  68   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
  69
  70   // Correction Terms for effective omegaTau; obtained by a laser calibration run
  71   Double_t T1 = 0.9;
  72   Double_t T2 = 1.5;
  73
  74   AliMagF mag("mag","mag");
  75
  76   AliTPCExBTwist twist;
  77   twist.SetXTwist(0.001);
  78
  79   AliTPCGGVoltError GGerror;
  80   GGerror.SetDeltaVGGA(50.);
  81   GGerror.SetDeltaVGGC(50.);
  82   GGerror.InitGGVoltErrorDistortion();
  83
  84   AliTPCExBBShape exb;
  85   exb.SetBField(&mag);
  86
  87   TObjArray cs;
  88   cs.Add(&twist);
  89   cs.Add(&GGerror);
  90   cs.Add(&exb);
  91
  92   AliTPCComposedCorrection cc;
  93   cc.SetCorrections(&cs);
  94   cc.SetOmegaTauT1T2(wt,T1,T2);
  95   //cc.SetMode(1);
  96
  97   cc.Print("DA"); // Print used correction classes
  98
  99   TCanvas *c=new TCanvas;  // Plots
 100   c->Divide(2,2);
 101   c->cd(1);twist.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");
 102   c->cd(2);GGerror.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");
 103   c->cd(3);exb.CreateHistoDRPhiinZR(1.)->Draw("surf2");
 104   c->cd(4);cc.CreateHistoDRPhiinZR(1.)->Draw("surf2");
 105 }
 106
 107
 108
 109 void MakeDistortionMap(){
 110   //
 111   // make distortiona map example for specific transformation
 112   //
 113   Int_t run=0;
 114   TTreeSRedirector *pcstream =  new TTreeSRedirector("distort.root");
 115   Double_t vdrift = 2.6; // [cm/us]   // to be updated: per second (ideally)
 116   Double_t bzField = -0.5; // [Tesla] // to be updated: per run
 117   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 118   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 119
 120   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 121   Double_t T1 = 0.9;
 122   Double_t T2 = 1.5;
 123   AliMagF* magF= new AliMagF("Maps","Maps", bzField/0.5, 1.);
 124   TGeoGlobalMagField::Instance()->SetField(magF);
 125
 126   AliCDBManager::Instance()->SetDefaultStorage("local:////$ALICE_ROOT/OCDB");
 127   AliCDBManager::Instance()->SetRun(run);
 128   AliGeomManager::LoadGeometry();
 129   AliGeomManager::ApplyAlignObjsFromCDB("GRP ITS TPC");
 130
 131
 132   AliTPCExBTwist twist;
 133   twist.SetName("twistX001");
 134   twist.SetXTwist(0.001);
 135
 136   AliTPCGGVoltError GGerror;
 137   GGerror.SetDeltaVGGA(50.);
 138   GGerror.SetDeltaVGGC(50.);
 139   GGerror.InitGGVoltErrorDistortion();
 140   GGerror.SetName("ggoffsetA50C50");
 141   AliTPCExBBShape exb;
 142   exb.SetBField(magF);
 143   exb.SetName("ExB");
 144   TObjArray cs;
 145   cs.Add(&twist);
 146   cs.Add(&GGerror);
 147   cs.Add(&exb);
 148
 149   AliTPCComposedCorrection *cc = new AliTPCComposedCorrection;
 150   cc->SetCorrections(&cs);
 151   cc->SetOmegaTauT1T2(wt,T1,T2);
 152   //cc->SetMode(1);
 153   cc->SetName("composed");
 154   Double_t par[5]={0,0,0,0,0};
 155   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 156   for (Double_t theta=-1.; theta<1; theta+=0.1){
 157     for (Double_t alpha=-3.14; alpha<3.14; alpha+=0.2){
 158       par[0]=0;
 159       par[1]=theta*85;
 160       par[2]=0;
 161       par[3]=theta;
 162       par[4]=(gRandom->Rndm()-0.5);
 163       AliExternalTrackParam trackIn(85,alpha,par,cov);
 164       //FitDistortedTrack(&trackIn, cc, 85, -1,pcstream);
 165       twist.FitDistortedTrack(trackIn, 85, -1,pcstream);
 166       GGerror.FitDistortedTrack(trackIn, 85, -1,pcstream);
 167       exb.FitDistortedTrack(trackIn, 85, -1,pcstream);
 168       cc->FitDistortedTrack(trackIn, 85, -1,pcstream);
 169     }
 170   }
 171   delete pcstream;
 172 }
 173
 174 void DrawDistortionMap(){
 175   //
 176   // Example -drawing of distortion maps
 177   //
 178   TFile f("distort.root");
 179   TTree *   fitDistorttwistX001= (TTree*)f.Get("fitDistorttwistX001");
 180   TTree *   fitDistortggoffsetA50C50=(TTree*)f.Get("fitDistortggoffsetA50C50");
 181   TTree *   fitDistortExB=(TTree*)f.Get("fitDistortExB");
 182   TTree *   fitDistortcomposed=(TTree*)f.Get("fitDistortcomposed");
 183   fitDistortExB->SetMarkerColor(1);
 184   fitDistorttwistX001->SetMarkerColor(2);
 185   fitDistortggoffsetA50C50->SetMarkerColor(4);
 186   fitDistortExB->SetMarkerStyle(20);
 187   fitDistorttwistX001->SetMarkerStyle(21);
 188   fitDistortggoffsetA50C50->SetMarkerStyle(22);
 189   //
 190   // example draw delta local y/r-phi as function of fi
 191   Int_t entries=0;
 192
 193   entries = fitDistortExB->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
 194   TGraph * grY0= new TGraph(entries, fitDistortExB->GetV2(),fitDistortExB->GetV1());
 195   entries = fitDistortggoffsetA50C50->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
 196   TGraph * grY1= new TGraph(entries, fitDistortggoffsetA50C50->GetV2(),fitDistortggoffsetA50C50->GetV1());
 197   entries = fitDistorttwistX001->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
 198   TGraph * grY2= new TGraph(entries, fitDistorttwistX001->GetV2(),fitDistorttwistX001->GetV1());
 199   entries = fitDistortcomposed->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
 200   TGraph * grY3= new TGraph(entries, fitDistortcomposed->GetV2(),fitDistortcomposed->GetV1());
 201   grY0->SetMinimum(-0.4);
 202   grY0->SetMaximum(0.4);
 203   grY0->SetMarkerStyle(20), grY0->SetMarkerColor(1);
 204   grY1->SetMarkerStyle(21), grY1->SetMarkerColor(2);
 205   grY2->SetMarkerStyle(22), grY2->SetMarkerColor(4);
 206   grY3->SetMarkerStyle(24), grY3->SetMarkerColor(5);
 207   grY0->Draw("ap");
 208   grY1->Draw("p");
 209   grY2->Draw("p");
 210   grY3->Draw("p");
 211 }
 212
 213
 214 void TestVertex(){
 215   //
 216   //
 217   //  .x ConfigCalibTrain.C(120829)
 218
 219   AliTPCComposedCorrection * corrC = ( AliTPCComposedCorrection *)AliTPCcalibDB::Instance()->GetTPCComposedCorrection(0.5);
 220   AliTPCCorrection * corrT = (AliTPCCorrection *)corrC->GetCorrections()->FindObject("exb_twist");
 221   AliTPCCorrection * corrS = (AliTPCCorrection *)corrC->GetCorrections()->FindObject("ExB");
 222
 223   TTreeSRedirector *pcstream = new TTreeSRedirector("vertexDistort.root");
 224   TTreeSRedirector *pcstreamS = new TTreeSRedirector("vertexDistortS.root");
 225   TTreeSRedirector *pcstreamT = new TTreeSRedirector("vertexDistortT.root");
 226   Double_t orgVertex[3] ;
 227   AliESDVertex aV,cV,aVO,cVO;
 228   for (Int_t iv=0; iv<100; iv++){
 229     printf("%d\n",iv);
 230     orgVertex[0]=gRandom->Gaus()*0.01;
 231     orgVertex[1]=gRandom->Gaus()*0.01;
 232     orgVertex[2]=gRandom->Gaus()*3;
 233     corrC->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstream);
 234     corrS->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstreamS);
 235     corrT->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstreamT);
 236   }
 237   delete pcstream;
 238   delete pcstreamS;
 239   delete pcstreamT;
 240   TFile fC("vertexDistortC.root");
 241   TFile fT("vertexDistortT.root");
 242   TFile fS("vertexDistortS.root");
 243   TTree * treeT=(TTree*)fT.Get("vertex"); // twist distrotions
 244   TTree * treeS=(TTree*)fS.Get("vertex"); // shape distortions
 245   TCanvas *canvasD=new TCanvas("canvasD","canvasD");
 246   canvasD->Divide(2,2);
 247   canvasD->cd(1);
 248   treeT->SetLineColor(2);
 249   treeT->Draw("av.fPosition[0]-x>>hisATX(100,-0.2,0.2)","","");
 250   treeT->SetLineColor(4);
 251   treeT->Draw("cv.fPosition[0]-x>>hisCTX(100,-0.2,0.2)","","same");
 252   canvasD->cd(2);
 253   treeT->SetLineColor(2);
 254   treeT->Draw("av.fPosition[1]-y>>hisATY(100,-0.2,0.2)","","");
 255   treeT->SetLineColor(4);
 256   treeT->Draw("cv.fPosition[1]-y>>hisCTY(100,-0.2,0.2)","","same");
 257   //
 258   canvasD->cd(3);
 259   treeS->SetLineColor(2);
 260   treeS->Draw("av.fPosition[0]-x>>hisASX(100,-0.2,0.2)","","");
 261   treeS->SetLineColor(4);
 262   treeS->Draw("cv.fPosition[0]-x>>hisCSX(100,-0.2,0.2)","","same");
 263   canvasD->cd(4);
 264   treeS->SetLineColor(2);
 265   treeS->Draw("av.fPosition[1]-y>>hisASY(100,-0.2,0.2)","","");
 266   treeS->SetLineColor(4);
 267   treeS->Draw("cv.fPosition[1]-y>>hisCSY(100,-0.2,0.2)","","same");
 268   canvasD->SaveAs("vertexShift.ps");
 269
 270 }