TPC/TPCupgrade/macros/SCdistortionsForTDRfinal.C

   1 void SCdistortionsForTDRfinalAll(){
   2
   3   SCdistortionsForTDRfinalsingle(1,1,0);
   4   SCdistortionsForTDRfinalsingle(1,1,1);
   5   SCdistortionsForTDRfinalsingle(1,1,2);
   6   SCdistortionsForTDRfinalsingle(1,1,3);
   7   SCdistortionsForTDRfinalsingle(1,1,4);
   8
   9 }
  10
  11 void SCdistortionsForTDRfinal(Double_t radiusScale=1.5, Int_t epsScale=1,Int_t iOmegaTau = 1){
  12   //
  13   // do for one file (given by directory) the correction (specify the gas = iOmegaTau)
  14   // use the integrate along drift line option
  15   //
  16   // 1. Initialzation form space charge maps
  17   //
  18   AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
  19   const Double_t fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
  20
  21   // omega tau parameters and TF1
  22   const Int_t nEps = 2;
  23   Int_t eps[nEps] = {20,10};
  24   Int_t col[nEps] = {kBlack,kBlack};
  25
  26   const Int_t nOmegaTau = 5;
  27   Double_t omegaTau[nOmegaTau] = {0.34,0.32,0.43,1.77,1.84};
  28   TString tGas[nOmegaTau] = {"NeCO2","NeCO2_2","ArCO2","NeCF4","NeCF4_2"}; // CF4 is the same as CO2 here, but different omegaTau
  29   TString sGas[nOmegaTau] = {"Ne-CO_{2} (90-10)","Ne-CO_{2}-N_{2} (90-10-5)","Ar-CO_{2} (90-10)","Ne-CF_{4} (90-10)","Ne-CF_{4} (80-20)"};
  30   TF2 * fdiffR[nEps];
  31   TF2 * fdiffPhiR[nEps];
  32   TH2F * hdiffR[nEps];
  33   TH2F * hdiffPhiR[nEps];
  34   TF2 * fdiffIntR[nEps];
  35   TF2 * fdiffIntPhiR[nEps];
  36   TF2 * fdiffIntZ[nEps];
  37   TH2F * hdiffIntR[nEps];
  38   TH2F * hdiffIntPhiR[nEps];
  39   TH2F * hdiffIntZ[nEps];
  40   TH2F * hMap[nEps];
  41   TH2F * hDistRMap[nEps];
  42   TH2F * hDistRPMap[nEps];
  43
  44   //use always the integrate option here
  45   Double_t integrateStep = 1.;
  46
  47   TCanvas *cMap = new TCanvas("cMap","cMap",1200,500);
  48   cMap->Divide(2,1);
  49
  50   TCanvas *cDistRMap = new TCanvas("cDistRMap","cDistRMap",1200,500);
  51   cDistRMap->Divide(2,1);
  52
  53   TCanvas *cDistRPMap = new TCanvas("cDistRPMap","cDistRPMap",1200,500);
  54   cDistRPMap->Divide(2,1);
  55
  56   TCanvas *cDistRNonIntMap = new TCanvas("cDistRNonIntMap","cDistRNonIntMap",1200,500);
  57   cDistRNonIntMap->Divide(2,1);
  58
  59   TCanvas *cDistRPNonIntMap = new TCanvas("cDistRPNonIntMap","cDistRPNonIntMap",1200,500);
  60   cDistRPNonIntMap->Divide(2,1);
  61
  62   TCanvas *cDistRIntMap = new TCanvas("cDistRIntMap","cDistRIntMap",1200,500);
  63   cDistRIntMap->Divide(2,1);
  64
  65   TCanvas *cDistRPIntMap = new TCanvas("cDistRPIntMap","cDistRPIntMap",1200,500);
  66   cDistRPIntMap->Divide(2,1);
  67
  68   TCanvas *cDistZIntMap = new TCanvas("cDistZIntMap","cDistZIntMap",1200,500);
  69   cDistZIntMap->Divide(2,1);
  70
  71
  72
  73   TString outfilename = Form("SCdistortions_%s_50kHz_radiusScaling%.0f_epsScaling%d",tGas[iOmegaTau].Data(),radiusScale,epsScale);
  74   if(radiusScale>1.1 && radiusScale < 1.9) outfilename = Form("SCdistortions_%s_50kHz_radiusScaling%.1f_epsScaling%d",tGas[iOmegaTau].Data(),radiusScale,epsScale);
  75
  76   //loop over epsilons
  77   for(Int_t iEps = 0; iEps < nEps; ++iEps){
  78
  79     cMap->cd(iEps+1);
  80
  81   // select gas
  82   // 0 = NeCO2
  83   // 1 = NeCO2N2
  84   // 2 = ArCO2
  85   // 3 = NeCF4
  86   // 4 = ArCF4
  87      if(radiusScale>1.1 && radiusScale < 1.9){
  88       cout<<"Open file = "<<Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.1f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale)<<endl;
  89       spaceCharge->SetSCDataFileName(Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.1f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale));
  90     }
  91     else{
  92       cout<<"Open file = "<<Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.0f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale)<<endl;
  93       spaceCharge->SetSCDataFileName(Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.0f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale));
  94     }
  95
  96   // select omegaTau value
  97   if(iOmegaTau ==  1){
  98     spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],1.00,0.99); // Ne CO2 N2 (90-10-5)
  99   }
 100   if(iOmegaTau ==  2){
 101     spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.99,1.03); // Ar CO2 (90-10)
 102   }
 103   else if(iOmegaTau == 3){
 104     spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.41,0.70); // Ne CF4 (90-10)
 105   }
 106   else if(iOmegaTau == 4){
 107     spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.41,0.70); // Ne CF4 (80-20) (not in table use same as for other)
 108   }
 109   else{
 110     spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],1,1.01); // Ne CO2 (90-10)
 111   }
 112   //
 113   //
 114   // init and add to corrections
 115   spaceCharge->InitSpaceCharge3DDistortion();
 116   spaceCharge->AddVisualCorrection(spaceCharge,1);
 117
 118   // draw the map
 119   hMap[iEps] = (TH2F*)spaceCharge->CreateHistoSCinZR(0.);
 120   hMap[iEps]->Scale(fgke0/1e6*1e15); // C/m^3/e0 --> fC/cm^3
 121   hMap[iEps]->GetXaxis()->SetTitle("z (cm)");
 122   hMap[iEps]->GetYaxis()->SetTitle("r (cm)");
 123   hMap[iEps]->GetZaxis()->SetTitle("");
 124   //hMap[iEps]->GetZaxis()->SetTitle("#rho_{SC} (fC/cm^{3})");
 125   hMap[iEps]->SetTitleSize(0.05,"XYZ");
 126   //hMap[iEps]->SetTitleOffset(1.5,"XY");
 127   //hMap[iEps]->SetTitleOffset(0.9,"Z");
 128   hMap[iEps]->SetTitle(Form("#rho_{SC} (fC/cm^{3}) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 129   hMap[iEps]->DrawCopy("colz");
 130
 131   // draw the distortion maps
 132   cDistRMap->cd(iEps+1);
 133   hDistRMap[iEps] = (TH2F*)spaceCharge->CreateHistoDRinZR(0.);
 134   hDistRMap[iEps]->SetMaximum(25);
 135   hDistRMap[iEps]->SetMinimum(-12);
 136   hDistRMap[iEps]->GetXaxis()->SetTitle("z (cm)");
 137   hDistRMap[iEps]->GetYaxis()->SetTitle("r (cm)");
 138   hDistRMap[iEps]->GetZaxis()->SetTitle("");
 139   hDistRMap[iEps]->SetTitleSize(0.05,"XYZ");
 140   //hDistRMap[iEps]->SetTitleOffset(1.5,"XY");
 141   //hDistRMap[iEps]->SetTitleOffset(0.9,"Z");
 142   hDistRMap[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 143   hDistRMap[iEps]->DrawCopy("colz");
 144
 145   cDistRPMap->cd(iEps+1);
 146   hDistRPMap[iEps] = (TH2F*)spaceCharge->CreateHistoDRinZR(0.);
 147   hDistRPMap[iEps]->SetMaximum(25);
 148   hDistRPMap[iEps]->SetMinimum(-12);
 149   hDistRPMap[iEps]->GetXaxis()->SetTitle("z (cm)");
 150   hDistRPMap[iEps]->GetYaxis()->SetTitle("r (cm)");
 151   hDistRPMap[iEps]->GetZaxis()->SetTitle("");
 152   hDistRPMap[iEps]->SetTitleSize(0.05,"XYZ");
 153   //hDistRPMap[iEps]->SetTitleOffset(1.5,"XY");
 154   //hDistRPMap[iEps]->SetTitleOffset(0.9,"Z");
 155   hDistRPMap[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 156   hDistRPMap[iEps]->DrawCopy("colz");
 157
 158
 159   //
 160   // 2. get TF2 with differences
 161   //
 162   // get corrections (at y = 0) for visual correction 1 (last argument)
 163   // 0 = dR
 164   // 1 = dPhiR
 165   // 2 = dZ?
 166   fdiffR[iEps]       = new TF2(Form("fdiffR%d",iEps), Form("AliTPCCorrection::GetDistXYZ(y,0,x,0,1)"),-250,250,85,250);
 167   fdiffPhiR[iEps]    = new TF2(Form("fdiffPhiR%d",iEps), Form("AliTPCCorrection::GetDistXYZ(y,0,x,1,1)"),-250,250,85,250);
 168   hdiffR[iEps] = (TH2F*)fdiffR[iEps]->GetHistogram();
 169   hdiffPhiR[iEps] = (TH2F*)fdiffPhiR[iEps]->GetHistogram();
 170
 171   hdiffR[iEps]->SetName(fdiffR[iEps]->GetName());
 172   hdiffPhiR[iEps]->SetName(fdiffPhiR[iEps]->GetName());
 173
 174   fdiffIntR[iEps]       = new TF2(Form("fdiffIntR%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,0,1,%f)",integrateStep),-250,250,85,250);
 175   fdiffIntPhiR[iEps]    = new TF2(Form("fdiffIntPhiR%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,1,1,%f)",integrateStep),-250,250,85,250);
 176   fdiffIntZ[iEps]    = new TF2(Form("fdiffIntZ%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,2,1,%f)",integrateStep),-250,250,85,250);
 177   hdiffIntR[iEps] = (TH2F*)fdiffIntR[iEps]->GetHistogram();
 178   hdiffIntPhiR[iEps] = (TH2F*)fdiffIntPhiR[iEps]->GetHistogram();
 179   hdiffIntZ[iEps] = (TH2F*)fdiffIntZ[iEps]->GetHistogram();
 180
 181   hdiffIntR[iEps]->SetName(fdiffIntR[iEps]->GetName());
 182   hdiffIntPhiR[iEps]->SetName(fdiffIntPhiR[iEps]->GetName());
 183   hdiffIntZ[iEps]->SetName(fdiffIntZ[iEps]->GetName());
 184
 185
 186
 187   //
 188   // 3. Plot and store TH1Fs
 189   //
 190
 191   // draw the distortion maps
 192   cDistRNonIntMap->cd(iEps+1);
 193   hdiffR[iEps]->SetMaximum(25);
 194   hdiffR[iEps]->SetMinimum(-12);
 195   hdiffR[iEps]->GetXaxis()->SetTitle("z (cm)");
 196   hdiffR[iEps]->GetYaxis()->SetTitle("r (cm)");
 197   //hdiffR[iEps]->GetZaxis()->SetTitle("dr (cm)");
 198   hdiffR[iEps]->SetTitleSize(0.05,"XYZ");
 199   //hdiffR[iEps]->SetTitleOffset(1.5,"XY");
 200   //hdiffR[iEps]->SetTitleOffset(0.9,"Z");
 201   hdiffR[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 202   hdiffR[iEps]->DrawCopy("colz");
 203
 204
 205   cDistRPNonIntMap->cd(iEps+1);
 206   hdiffPhiR[iEps]->SetMaximum(25);
 207   hdiffPhiR[iEps]->SetMinimum(-12);
 208   hdiffPhiR[iEps]->GetXaxis()->SetTitle("z (cm)");
 209   hdiffPhiR[iEps]->GetYaxis()->SetTitle("r (cm)");
 210   hdiffPhiR[iEps]->SetTitleSize(0.05,"XYZ");
 211   //hdiffPhiR[iEps]->SetTitleOffset(1.5,"XY");
 212   //hdiffPhiR[iEps]->SetTitleOffset(0.9,"Z");
 213   hdiffPhiR[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 214   hdiffPhiR[iEps]->DrawCopy("colz");
 215
 216   cDistRIntMap->cd(iEps+1);
 217   hdiffIntR[iEps]->SetMaximum(25);
 218   hdiffIntR[iEps]->SetMinimum(-12);
 219   hdiffIntR[iEps]->GetXaxis()->SetTitle("z (cm)");
 220   hdiffIntR[iEps]->GetYaxis()->SetTitle("r (cm)");
 221   //hdiffIntR[iEps]->GetZaxis()->SetTitle("dr (cm)");
 222   hdiffIntR[iEps]->SetTitleSize(0.05,"XYZ");
 223   //hdiffIntR[iEps]->SetTitleOffset(1.5,"XY");
 224   //hdiffIntR[iEps]->SetTitleOffset(0.9,"Z");
 225   hdiffIntR[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 226   hdiffIntR[iEps]->DrawCopy("colz");
 227
 228
 229   cDistRPIntMap->cd(iEps+1);
 230   hdiffIntPhiR[iEps]->SetMaximum(5);
 231   hdiffIntPhiR[iEps]->SetMinimum(-9);
 232   hdiffIntPhiR[iEps]->GetXaxis()->SetTitle("z (cm)");
 233   hdiffIntPhiR[iEps]->GetYaxis()->SetTitle("r (cm)");
 234   hdiffIntPhiR[iEps]->SetTitleSize(0.05,"XYZ");
 235   //hdiffIntPhiR[iEps]->SetTitleOffset(1.5,"XY");
 236   //hdiffIntPhiR[iEps]->SetTitleOffset(0.9,"Z");
 237   hdiffIntPhiR[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 238   hdiffIntPhiR[iEps]->DrawCopy("colz");
 239
 240   cDistZIntMap->cd(iEps+1);
 241   hdiffIntZ[iEps]->SetMaximum(5);
 242   hdiffIntZ[iEps]->SetMinimum(-5);
 243   hdiffIntZ[iEps]->GetXaxis()->SetTitle("z (cm)");
 244   hdiffIntZ[iEps]->GetYaxis()->SetTitle("r (cm)");
 245   hdiffIntZ[iEps]->SetTitleSize(0.05,"XYZ");
 246   //hdiffIntZ[iEps]->SetTitleOffset(1.5,"XY");
 247   //hdiffIntZ[iEps]->SetTitleOffset(0.9,"Z");
 248   hdiffIntZ[iEps]->SetTitle(Form("dz (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
 249   hdiffIntZ[iEps]->DrawCopy("colz");
 250   }
 251
 252   // just control histograms (not saved)
 253   // cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.eps",outfilename.Data()));
 254   // cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.png",outfilename.Data()));
 255   // cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.pdf",outfilename.Data()));
 256   // cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.eps",outfilename.Data()));
 257   // cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.png",outfilename.Data()));
 258   // cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.pdf",outfilename.Data()));
 259
 260   // histograms for TDR (saved)
 261   cMap->SaveAs(Form("%s_TDR_SpaceCharge.eps",outfilename.Data()));
 262   cMap->SaveAs(Form("%s_TDR_SpaceCharge.png",outfilename.Data()));
 263   cMap->SaveAs(Form("%s_TDR_SpaceCharge.pdf",outfilename.Data()));
 264   cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.eps",outfilename.Data()));
 265   cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.png",outfilename.Data()));
 266   cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.pdf",outfilename.Data()));
 267   cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.eps",outfilename.Data()));
 268   cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.png",outfilename.Data()));
 269   cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.pdf",outfilename.Data()));
 270   cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.eps",outfilename.Data()));
 271   cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.png",outfilename.Data()));
 272   cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.pdf",outfilename.Data()));
 273 }
 274
 275
 276 //____________________________________________________________//
 277 void setupLegend(TLegend *currentLegend=0,float currentTextSize=0.07){
 278   currentLegend->SetTextFont(42);
 279   currentLegend->SetBorderSize(0);
 280   currentLegend->SetFillStyle(0);
 281   currentLegend->SetFillColor(0);
 282   currentLegend->SetMargin(0.25);
 283   currentLegend->SetTextSize(currentTextSize);
 284   currentLegend->SetEntrySeparation(0.5);
 285   return;
 286 }