TPC/CalibMacros/AliTPCCorrectionDemo.C

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