]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMagF.cxx
Implementation of a method (RemoveGraphDuplicates) to take away
[u/mrichter/AliRoot.git] / STEER / AliMagF.cxx
index fa20c79a01fd5a47b2bfc04e2ef6377002dbd0b0..f7960fece85c5cf9392810c8ab378abb62a3721d 100644 (file)
@@ -17,6 +17,7 @@
 #include <TClass.h>
 #include <TFile.h>
 #include <TSystem.h>
+#include <TPRegexp.h>
 
 #include "AliMagF.h"
 #include "AliMagWrapCheb.h"
 ClassImp(AliMagF)
 
 const Double_t AliMagF::fgkSol2DipZ    =  -700.;  
+const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
 
+/*
+ Explanation for polarity conventions: these are the mapping between the
+ current signs and main field components in L3 (Bz) and Dipole (Bx) (in Alice frame)
+ 1) kConvMap2005: used for the field mapping in 2005
+ positive L3  current -> negative Bz
+ positive Dip current -> positive Bx 
+ 2) kConvMapDCS2008: defined by the microswitches/cabling of power converters as of 2008 - 1st half 2009
+ positive L3  current -> positive Bz
+ positive Dip current -> positive Bx
+ 3) kConvLHC : defined by LHC
+ positive L3  current -> positive Bz
+ positive Dip current -> negative Bx
+ Note: only "negative Bz(L3) with postive Bx(Dipole)" and its inverse was mapped in 2005. Hence 
+ the GRP Manager will reject the runs with the current combinations (in the convention defined by the
+ static Int_t AliMagF::GetPolarityConvention()) which do not lead to such field polarities.
+*/
 //_______________________________________________________________________
 AliMagF::AliMagF():
   TVirtualMagField(),
@@ -54,10 +73,8 @@ AliMagF::AliMagF():
 }
 
 //_______________________________________________________________________
-AliMagF::AliMagF(const char *name, const char* title, Int_t integ, 
-                Double_t factorSol, Double_t factorDip, 
-                Double_t fmax, BMap_t maptype, const char* path,
-                BeamType_t bt, Double_t be):
+AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip, 
+                BMap_t maptype, BeamType_t bt, Double_t be,Int_t integ, Double_t fmax, const char* path):
   TVirtualMagField(name),
   fMeasuredMap(0),
   fMapType(maptype),
@@ -90,28 +107,32 @@ AliMagF::AliMagF(const char *name, const char* title, Int_t integ,
   }
   if (fInteg == 0) fPrecInteg = 0;
   //
+  if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
+    if      (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
+    else if (fBeamType == kBeamTypeAA) fBeamEnergy = 2750;  // max PbPb energy
+    AliInfo("Maximim possible beam energy for requested beam is assumed");
+  } 
   const char* parname = 0;
   //  
-  if (fMapType == k2kG) {
-    fSolenoid = 2.;
-    parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
-  } else if (fMapType == k5kG) {
-    fSolenoid = 5.;
-    parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
-  } else if (fMapType == k5kGUniform) {
-    fSolenoid = 5.;
-    parname = "Sol30_Dip6_Uniform";
-  } else {
-    AliFatal(Form("Unknown field identifier %d is requested\n",fMapType)); 
-  }
+  if      (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
+  else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
+  else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
+  else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
   //
   SetDataFileName(path);
   SetParamName(parname);
   //
-  SetFactorSol(factorSol);
-  SetFactorDip(factorDip);
   LoadParameterization();
   InitMachineField(fBeamType,fBeamEnergy);
+  double xyz[3]={0.,0.,0.};
+  fSolenoid = GetBz(xyz);
+  SetFactorSol(factorSol);
+  SetFactorDip(factorDip);
+  AliInfo(Form("Alice   B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
+              factorSol,(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
+              fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" |Constant Field!":""));
+  AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
+              bt==kBeamTypeAA ? "A-A":(bt==kBeamTypepp ? "p-p":"OFF"),be,fQuadGradient,fDipoleField));
 }
 
 //_______________________________________________________________________
@@ -221,7 +242,7 @@ AliMagF& AliMagF::operator=(const AliMagF& src)
 //_______________________________________________________________________
 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
 {
-  if (btype==kNoBeamField || benergy<1.) {
+  if (btype==kNoBeamField) {
     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
     return;
   }
@@ -250,7 +271,7 @@ void AliMagF::MachineField(const Double_t *x, Double_t *b) const
   const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; 
   //  
   const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
-  const Double_t kTripQ2CZ = 3408., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
+  const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
   const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
   const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
   //
@@ -355,3 +376,144 @@ void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
     for (int i=3;i--;) b[i] *= fFactorSol;
   }
 }
+
+//_______________________________________________________________________
+void AliMagF::SetFactorSol(Float_t fc)
+{
+  // set the sign/scale of the current in the L3 according to fgkPolarityConvention
+  switch (fgkPolarityConvention) {
+  case kConvDCS2008: fFactorSol = -fc; break;
+  case kConvLHC    : fFactorSol = -fc; break;
+  default          : fFactorSol =  fc; break;  // case kConvMap2005: fFactorSol =  fc; break;
+  }
+}
+
+//_______________________________________________________________________
+void AliMagF::SetFactorDip(Float_t fc)
+{
+  // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
+  switch (fgkPolarityConvention) {
+  case kConvDCS2008: fFactorDip =  fc; break;
+  case kConvLHC    : fFactorDip = -fc; break;
+  default          : fFactorDip =  fc; break;  // case kConvMap2005: fFactorDip =  fc; break;
+  }
+}
+
+//_______________________________________________________________________
+Double_t AliMagF::GetFactorSol() const
+{
+  // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
+  switch (fgkPolarityConvention) {
+  case kConvDCS2008: return -fFactorSol;
+  case kConvLHC    : return -fFactorSol;
+  default          : return  fFactorSol;       //  case kConvMap2005: return  fFactorSol;
+  }
+}
+
+//_______________________________________________________________________
+Double_t AliMagF::GetFactorDip() const
+{
+  // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
+  switch (fgkPolarityConvention) {
+  case kConvDCS2008: return  fFactorDip;
+  case kConvLHC    : return -fFactorDip;
+  default          : return  fFactorDip;       //  case kConvMap2005: return  fFactorDip;
+  }
+}
+
+//_____________________________________________________________________________
+AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
+                                Float_t beamenergy, const Char_t *beamtype, const Char_t *path) 
+{
+  //------------------------------------------------
+  // The magnetic field map, defined externally...
+  // L3 current 30000 A  -> 0.5 T
+  // L3 current 12000 A  -> 0.2 T
+  // dipole current 6000 A
+  // The polarities must match the convention (LHC or DCS2008) 
+  // unless the special uniform map was used for MC
+  //------------------------------------------------
+  const Float_t l3NominalCurrent1=30000.; // (A)
+  const Float_t l3NominalCurrent2=12000.; // (A)
+  const Float_t diNominalCurrent =6000. ; // (A)
+
+  const Float_t tolerance=0.03; // relative current tolerance
+  const Float_t zero=77.;       // "zero" current (A)
+  //
+  BMap_t map;
+  double sclL3,sclDip;
+  //
+  Float_t l3Pol = l3Cur > 0 ? 1:-1;
+  Float_t diPol = diCur > 0 ? 1:-1;
+  l3Cur = TMath::Abs(l3Cur);
+  diCur = TMath::Abs(diCur);
+  //
+  if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
+    if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
+    else {
+      AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
+      return 0;
+    }
+  }
+  //
+  if (uniform) { 
+    // special treatment of special MC with uniform mag field (normalized to 0.5 T)
+    // no check for scaling/polarities are done
+    map   = k5kGUniform;
+    sclL3 = l3Cur/l3NominalCurrent1; 
+  }
+  else {
+    if      (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map  = k5kG;
+    else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map  = k2kG;
+    else if (l3Cur <= zero && diCur<=zero)   { sclL3=0; sclDip=0; map  = k5kGUniform;}
+    else {
+      AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
+      return 0;
+    }
+  }
+  //
+  if (sclDip!=0 && map!=k5kGUniform) {
+    if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) { 
+      AliErrorGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
+                                    l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
+      return 0;
+    }
+  }
+  //
+  if (l3Pol<0) sclL3  = -sclL3;
+  if (diPol<0) sclDip = -sclDip;
+  //
+  BeamType_t btype = kNoBeamField;
+  TString btypestr = beamtype;
+  btypestr.ToLower();
+  TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
+  TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1");
+  if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
+  else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
+  else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
+  char ttl[80];
+  sprintf(ttl,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
+         (int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"",
+         convention==kConvLHC ? "LHC":"DCS2008");
+  // LHC and DCS08 conventions have opposite dipole polarities
+  if ( GetPolarityConvention() != convention) sclDip = -sclDip;
+  //
+  return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
+  //
+}
+
+//_____________________________________________________________________________
+const char*  AliMagF::GetBeamTypeText() const
+{
+  const char *beamNA  = "No Beam";
+  const char *beamPP  = "p-p";
+  const char *beamPbPb= "Pb-Pb";
+  switch ( fBeamType ) {
+  case kBeamTypepp : return beamPP;
+  case kBeamTypeAA : return beamPbPb;
+  case kNoBeamField: 
+  default:           return beamNA;
+  }
+}
+