TPC/TPCupgrade/macros/spaceChargeSimulationTrack.C

   1 void CheckCommutator(const char * inputCharge="SpaceCharge.root"){
   2   //
   3   // 1.) create an space charge correction map - using as input 2x2d space charge distribution
   4   // file
   5   /*
   6     const char * inputCharge="SpaceCharge.root"
   7   */
   8   AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
   9   spaceCharge->SetSCDataFileName(inputCharge);
  10   spaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
  11   //spaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
  12   spaceCharge->InitSpaceCharge3DDistortion();
  13   spaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
  14   spaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
  15   spaceCharge->AddVisualCorrection(spaceCharge,1);
  16   //
  17   // 2.) Instantiate magnetic field
  18   //
  19   ocdb="local://$ALICE_ROOT/OCDB/";
  20   AliCDBManager::Instance()->SetDefaultStorage(ocdb);
  21   AliCDBManager::Instance()->SetRun(0);
  22   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  23   //
  24   // 3.)   OCDB correction  have to be initialized from somewhere
  25   //
  26   TFile * f  = TFile::Open("Correction.root");
  27   TObjArray  * array = AliCDBEntry->GetObject();
  28   AliTPCComposedCorrection * corrDefault = (AliTPCComposedCorrection *)array->At(0);
  29   corrDefault->Init();
  30   //
  31   // 4.) Create comutators
  32   //
  33   TObjArray * array_SC_Default  = new TObjArray(2);
  34   TObjArray * array_Default_SC  = new TObjArray(2);
  35   //
  36   array_SC_Default->AddAt(spaceCharge,0);
  37   array_SC_Default->AddAt(corrDefault,1);
  38   //
  39   array_Default_SC->AddAt(corrDefault,0);
  40   array_Default_SC->AddAt(spaceCharge,1);
  41    AliTPCComposedCorrection *corr_SC_Default = new AliTPCComposedCorrection(array_SC_Default,AliTPCComposedCorrection::kQueue);
  42   AliTPCComposedCorrection *corr_Default_SC = new AliTPCComposedCorrection(array_Default_SC,AliTPCComposedCorrection::kQueue);
  43   corr_SC_Default->AddVisualCorrection(corr_SC_Default,1);
  44   corr_Default_SC->AddVisualCorrection(corr_Default_SC,2);
  45   //
  46   // 5.) Use TF1, TF2 functionality to visualize comutators
  47   //
  48   //
  49   TCanvas * canvasNonComute = new TCanvas("SpacechargeDefault"," SpacechargeDefault",1200,700);
  50   canvasNonComute->Divide(2,1);
  51   canvasNonComute->cd(1);
  52
  53   canvasNonComute->cd(1);gPad->SetRightMargin(0.2);
  54   TF2 * fdiffXY0 = new TF2("fdiffXY0", "(AliTPCCorrection::GetCorrXYZ(x,y,10,0,1)-AliTPCCorrection::GetCorrXYZ(x,y,10,0,2))*(sqrt(x*x+y*y)>85&&sqrt(x*x+y*y)<245)",-250,250,-250,250);
  55   fdiffXY0->SetNpx(200);
  56   fdiffXY0->SetNpy(200);
  57   fdiffXY0->GetXaxis()->SetTitle("x (cm)");
  58   fdiffXY0->GetYaxis()->SetTitle("y (cm)");
  59   fdiffXY0->GetZaxis()->SetTitle("#delta R (cm)");
  60   fdiffXY0->Draw("colz");
  61   TPaveText *paveText = new TPaveText(-250,200,0.0,250,"");
  62   paveText->AddText("[SpaceCharge,Default]");
  63   paveText->Draw();
  64
  65
  66   canvasNonComute->cd(2);gPad->SetRightMargin(0.2);
  67   TF2 * fdiffXY1 = new TF2("fdiffXY0", "(AliTPCCorrection::GetCorrXYZ(x,y,10,1,1)-AliTPCCorrection::GetCorrXYZ(x,y,10,1,2))*(sqrt(x*x+y*y)>85&&sqrt(x*x+y*y)<245)",-250,250,-250,250);
  68   fdiffXY1->SetNpx(200);
  69   fdiffXY1->SetNpy(200);
  70   fdiffXY1->GetXaxis()->SetTitle("x (cm)");
  71   fdiffXY1->GetYaxis()->SetTitle("y (cm)");
  72   fdiffXY1->GetZaxis()->SetTitle("#delta R#Phi (cm)");
  73   fdiffXY1->Draw("colz");
  74   //
  75   //
  76   //
  77 }
  78
  79
  80
  81
  82 void DrawFuncionIntegralSpaceCharge(){
  83   //
  84   // Make nice plot and compare the r and rphi distortion using "Stefan standard implementation" integration
  85   // and using integration along trajectories
  86   //
  87   AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
  88   spaceCharge->SetSCDataFileName("SpaceCharge.root");
  89   spaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
  90   //spaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
  91   spaceCharge->InitSpaceCharge3DDistortion();
  92   spaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
  93   spaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
  94   ocdb="local://$ALICE_ROOT/OCDB/";
  95   AliCDBManager::Instance()->SetDefaultStorage(ocdb);
  96   AliCDBManager::Instance()->SetRun(0);
  97   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  98   spaceCharge->AddVisualCorrection(spaceCharge,1);
  99
 100   //
 101   //
 102   TCanvas *canvasIntegrate = new TCanvas("canvasIntegrate","canvasIntegrate",600,600);
 103   canvasIntegrate->Divide(1,2);
 104   //
 105   TF1 *fdistRStefan = new TF1("fdistRStefan","AliTPCCorrection::GetCorrXYZ(x,0,10,0,1)",85,250);
 106   TF1 *fdistRDriftS5 = new TF1("fdistRDriftS3","AliTPCCorrection::GetCorrXYZIntegrateZ(x,0,10,0,1,5)",85,250);
 107   TF1 *fdistRDriftS2 = new TF1("fdistRDriftS2","AliTPCCorrection::GetCorrXYZIntegrateZ(x,0,10,0,1,2)",85,250);
 108   TF1 *fdistRDriftS1 = new TF1("fdistRDriftS1","AliTPCCorrection::GetCorrXYZIntegrateZ(x,0,10,0,1,1)",85,250);
 109   TF1 *fdistRPhiStefan = new TF1("fdistRPhiStefan","AliTPCCorrection::GetCorrXYZ(x,0,10,1,1)",85,250);
 110   TF1 *fdistRPhiDriftS5 = new TF1("fdistRPhiDriftS3","AliTPCCorrection::GetCorrXYZIntegrateZ(x,0,10,1,1,5)",85,250);
 111   TF1 *fdistRPhiDriftS2 = new TF1("fdistRPhiDriftS2","AliTPCCorrection::GetCorrXYZIntegrateZ(x,0,10,1,1,2)",85,250);
 112   TF1 *fdistRPhiDriftS1 = new TF1("fdistRPhiDriftS1","AliTPCCorrection::GetCorrXYZIntegrateZ(x,0,10,1,1,1)",85,250);
 113   canvasIntegrate->cd(1);
 114   TLegend * legendR=new TLegend(0.4,0.6,0.9,0.9,"Space charge #Delta_{R}");
 115   fdistRStefan->GetXaxis()->SetTitle("R (cm)");
 116   fdistRStefan->GetYaxis()->SetTitle("#Delta_{R} (cm)");
 117   fdistRStefan->SetLineColor(2);
 118   //
 119   fdistRDriftS5->SetLineColor(1); fdistRDriftS5->SetLineStyle(2);
 120   fdistRDriftS2->SetLineColor(4); fdistRDriftS2->SetLineStyle(1);
 121   fdistRDriftS1->SetLineColor(1); fdistRDriftS1->SetLineStyle(3);
 122   fdistRStefan->Draw();
 123   fdistRDriftS5->Draw("same");
 124   fdistRDriftS2->Draw("same");
 125   fdistRDriftS1->Draw("same");
 126   legendR->AddEntry(fdistRStefan,"Stefan integration");
 127   legendR->AddEntry(fdistRDriftS5,"Drift lines (step=5cm);");
 128   legendR->AddEntry(fdistRDriftS2,"Drift lines (step=2cm);");
 129   legendR->AddEntry(fdistRDriftS1,"Drift lines (step=1cm);");
 130   legendR->Draw();
 131
 132   canvasIntegrate->cd(2);
 133   fdistRPhiStefan->GetXaxis()->SetTitle("RPhi (cm)");
 134   fdistRPhiStefan->GetYaxis()->SetTitle("#Delta_{RPhi} (cm)");
 135   fdistRPhiStefan->SetLineColor(2);
 136   //
 137   fdistRPhiDriftS5->SetLineColor(4); fdistRPhiDriftS5->SetLineStyle(2);
 138   fdistRPhiDriftS2->SetLineColor(4);
 139   fdistRPhiDriftS1->SetLineColor(4); fdistRPhiDriftS1->SetLineStyle(2);
 140   fdistRPhiStefan->Draw();
 141   fdistRPhiDriftS5->Draw("same");
 142   fdistRPhiDriftS2->Draw("same");
 143   fdistRPhiDriftS1->Draw("same");
 144
 145   TLegend * legendRPhi=new TLegend(0.4,0.6,0.9,0.9,"Space charge #Delta_{RPhi}");
 146   fdistRPhiStefan->GetXaxis()->SetTitle("RPhi (cm)");
 147   fdistRPhiStefan->GetYaxis()->SetTitle("#Delta_{RPhi} (cm)");
 148   fdistRPhiStefan->SetLineColor(2);
 149  legendRPhi->AddEntry(fdistRPhiStefan,"Stefan integration");
 150   legendRPhi->AddEntry(fdistRPhiDriftS5,"Drift lines (step=5cm);");
 151   legendRPhi->AddEntry(fdistRPhiDriftS2,"Drift lines (step=2cm);");
 152   legendRPhi->AddEntry(fdistRPhiDriftS1,"Drift lines (step=1cm);");
 153   legendRPhi->Draw();
 154   canvasIntegrate->SaveAs("canvasIntegrate.png");
 155 }
 156
 157 void DrawFuncions(){
 158   //
 159   // Make a default plot for the
 160   //
 161   AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
 162   spaceCharge->SetSCDataFileName("SpaceCharge.root");
 163   spaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
 164   //spaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
 165   spaceCharge->InitSpaceCharge3DDistortion();
 166   spaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
 167   spaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
 168   ocdb="local://$ALICE_ROOT/OCDB/";
 169   AliCDBManager::Instance()->SetDefaultStorage(ocdb);
 170   AliCDBManager::Instance()->SetRun(0);
 171   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 172   spaceCharge->AddVisualCorrection(spaceCharge,1);
 173
 174   TCanvas *canvasFun = new TCanvas("canvasFun","canvasFun",600,500);
 175   gStyle->SetOptTitle(1);
 176   TF1 * fdiffR = new TF1("fdiffR", "(AliTPCCorrection::GetCorrXYZ(x,0,10,0,1)-AliTPCCorrection::GetCorrXYZ(x+1,0,10,0,1))",80,245);
 177   fdiffR->SetNpx(1000);
 178   fdiffR->GetXaxis()->SetTitle("R (cm)");
 179   fdiffR->GetYaxis()->SetTitle("#Delta_{R}(R)-#Delta_{R}(R+1) (cm)");
 180   fdiffR->Draw();
 181   canvasFun->SaveAs("radialShrinking.png");
 182
 183   TF1 * fdiffSigmaR = new TF1("fdiffSigmaR", "(AliTPCCorrection::GetCorrXYZ(x,0,10,0,1)-AliTPCCorrection::GetCorrXYZ(x,0,10-7,0,1))",80,245);
 184   fdiffSigmaR->SetNpx(1000);
 185   fdiffSigmaR->GetXaxis()->SetTitle("R (cm)");
 186   fdiffSigmaR->GetYaxis()->SetTitle("#Delta_{R}(R,Z)-#Delta_{R}(R,Z+#sigma_{z}) (cm)");
 187   fdiffSigmaR->Draw();
 188   canvasFun->SaveAs("radialSigmaR.png");
 189   //
 190   TF1 * fdistortion = new TF1("fdiffSigmaR", "AliTPCCorrection::GetCorrXYZ(x,0,10,0,1)+x",80,245);
 191   fdistortion->SetNpx(1000);
 192   fdistortion->GetXaxis()->SetTitle("R (cm)");
 193   fdistortion->GetYaxis()->SetTitle("R' (cm)");
 194   fdistortion->Draw();
 195   canvasFun->SaveAs("radialDistortion.png");
 196 }
 197
 198
 199 void spaceChargeSimulationTrack(Double_t refX){
 200   //
 201   // 1. Initialzation form space charge maps
 202   //
 203   AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
 204   spaceCharge->SetSCDataFileName("SpaceCharge.root");
 205   spaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
 206   //spaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
 207   spaceCharge->InitSpaceCharge3DDistortion();
 208   spaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
 209   spaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
 210   ocdb="local://$ALICE_ROOT/OCDB/";
 211   AliCDBManager::Instance()->SetDefaultStorage(ocdb);
 212   AliCDBManager::Instance()->SetRun(0);
 213   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 214   spaceCharge->AddVisualCorrection(spaceCharge,1);
 215   TF1 * fdiffR = new TF1("fdiffR", "(AliTPCCorrection::GetCorrXYZ(x,0,10,0,1)-AliTPCCorrection::GetCorrXYZ(x+1,0,10,0,1))",80,245);
 216   //
 217   // 2. MC generation of the tracks
 218   //    Distort track and fit distroted track within AliTPCorrection::FitDistortedTrack(*t, refX, dir,  pcstream)
 219   //
 220   Double_t etaCuts=1;
 221   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 222   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
 223   fpt.SetParameters(7.24,0.120);
 224   fpt.SetNpx(10000);
 225   Int_t nTracks=10000;
 226   TTreeSRedirector  * pcstream = new TTreeSRedirector("trackDist.root","recreate");
 227
 228   for(Int_t nt=0; nt<nTracks; nt++){
 229     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
 230     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
 231     Double_t pt = 1/(gRandom->Rndm()*5+0.00001); // momentum for f1
 232     //   printf("phi %lf  eta %lf pt %lf\n",phi,eta,pt);
 233     Short_t sign=1;
 234     if(gRandom->Rndm() < 0.5){
 235       sign =1;
 236     }else{
 237       sign=-1;
 238     }
 239
 240     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
 241     Double_t pxyz[3];
 242     pxyz[0]=pt*TMath::Cos(phi);
 243     pxyz[1]=pt*TMath::Sin(phi);
 244     pxyz[2]=pt*TMath::Tan(theta);
 245     Double_t vertex[3]={0,0,0};
 246     Double_t cv[21]={0};
 247     AliExternalTrackParam *t= new AliExternalTrackParam(vertex, pxyz, cv, sign);
 248     //    Double_t refX=1.;
 249     Int_t dir=-1;
 250     AliExternalTrackParam *td =  spaceCharge->FitDistortedTrack(*t, refX, dir,  pcstream);
 251   }
 252   delete pcstream;
 253   //
 254   // Visulalize distortion of of tracks
 255   //
 256   TFile * f = TFile::Open("trackDist.root");
 257   TTree * tree = (TTree*)f->Get("fitDistortSpaceCharge3D");
 258   tree->SetMarkerStyle(25);
 259   tree->SetMarkerSize(0.4);
 260
 261   //
 262   // DCA distortion example
 263   //
 264   tree->Draw("track1.GetY()-track0.GetY():track0.GetTgl():abs(track0.fP[4])","track0.fP[4]>0","colz",10000);
 265   tree->Draw("track1.GetTgl()-track0.GetTgl():track0.GetSigned1Pt():track0.GetTgl()","abs(track0.GetTgl()-0.5)<0.4","colz",10000);
 266
 267   /*
 268     problems:
 269
 270    */
 271 }
 272
 273 /*
 274  */
 275
 276 void DrawDistortions(){
 277   //
 278   //
 279   //
 280   TFile * f = TFile::Open("trackDist.root");
 281   TTree * tree = (TTree*)f->Get("fitDistortSpaceCharge3D");
 282   tree->SetMarkerStyle(25);
 283   tree->SetMarkerSize(0.4);
 284
 285   //
 286   // DCA as function of theta and 1/pt
 287   //
 288   TCanvas *c1 = new TCanvas();
 289   tree->Draw("track1.GetY()-track0.GetY():track0.GetSigned1Pt():track0.GetTgl()","track0.GetTgl()>0","colz",10000);
 290   // Get the histogram and set axis titles
 291   TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
 292   htemp->GetXaxis()->SetTitle("1/p_{T}");
 293   htemp->GetYaxis()->SetTitle("#Delta DCA");
 294   htemp->GetZaxis()->SetTitle("tan(#lambda)");
 295   htemp->SetTitle("DCA difference");
 296   c1->Update();
 297   c1->SaveAs("DCAdifference.eps")
 298
 299   //
 300   tree->Draw("track1.GetSigned1Pt()-track0.GetSigned1Pt():track0.GetSigned1Pt():track0.GetTgl()","track0.GetTgl()>0","colz",10000);
 301   htemp = (TH2F*)gPad->GetPrimitive("htemp");
 302   htemp->GetXaxis()->SetTitle("1/p_{T}");
 303   htemp->GetYaxis()->SetTitle("#Delta 1/pT (1/GeV/c)");
 304   htemp->GetZaxis()->SetTitle("tan(#lambda)");
 305   htemp->SetTitle("1/pT difference");
 306   c1->Update();
 307   c1->SaveAs("1Ptdifference.eps")
 308   //
 309   //
 310   TCanvas *canvasEvent = new TCanvas("canvasEvent","canvasEvent",600,500);
 311   tree->SetMarkerSize(0.25);
 312   tree->SetMarkerColor(1);
 313   tree->Draw("point0.fY:point0.fX","point0.fX!=0","",40,0);
 314   tree->SetMarkerColor(2);
 315   tree->Draw("point1.fY:point1.fX","point1.fX!=0","same",40,0);
 316   htemp = (TH2F*)gPad->GetPrimitive("htemp");
 317   htemp->GetXaxis()->SetTitle("x (cm)");
 318   htemp->GetYaxis()->SetTitle("y (cm)");
 319   htemp->SetTitle("Event");
 320   c1->Update();
 321   c1->SaveAs("Event.eps")
 322
 323 }