]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMagF.cxx
New centrality files. (A. Toia)
[u/mrichter/AliRoot.git] / STEER / AliMagF.cxx
index c32309ee256209ab5f9617ce3688a2fd6403481a..e709b6ef71f3408e96346e2c1902be93ed86bcd2 100644 (file)
@@ -26,8 +26,7 @@
 ClassImp(AliMagF)
 
 const Double_t AliMagF::fgkSol2DipZ    =  -700.;  
-const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
-
+const UShort_t AliMagF::fgkPolarityConvention = AliMagF::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)
@@ -38,12 +37,31 @@ const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
  positive L3  current -> positive Bz
  positive Dip current -> positive Bx
  3) kConvLHC : defined by LHC
- positive L3  current -> negative Bz
+ 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.
+
+ ----------------------------------------------- 
+
+ Explanation on integrals in the TPC region
+ GetTPCInt(xyz,b) and GetTPCRatInt(xyz,b) give integrals from point (x,y,z) to point (x,y,0) 
+ (irrespectively of the z sign) of the following:
+ TPCInt:    b contains int{bx}, int{by}, int{bz}
+ TPCRatInt: b contains int{bx/bz}, int{by/bz}, int{(bx/bz)^2+(by/bz)^2}
+  
+ The same applies to integral in cylindrical coordinates:
+ GetTPCIntCyl(rphiz,b)
+ GetTPCIntRatCyl(rphiz,b)
+ They accept the R,Phi,Z coordinate (-pi<phi<pi) and return the field 
+ integrals in cyl. coordinates.
+
+ Thus, to compute the integral from arbitrary xy_z1 to xy_z2, one should take
+ b = b1-b2 with b1 and b2 coming from GetTPCInt(xy_z1,b1) and GetTPCInt(xy_z2,b2)
+
+ Note: the integrals are defined for the range -300<Z<300 and 0<R<300
 */
 //_______________________________________________________________________
 AliMagF::AliMagF():
@@ -73,10 +91,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),
@@ -109,6 +125,11 @@ 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) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
@@ -125,11 +146,7 @@ AliMagF::AliMagF(const char *name, const char* title, Int_t integ,
   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));
+  Print("a");
 }
 
 //_______________________________________________________________________
@@ -166,21 +183,18 @@ AliMagF::~AliMagF()
 Bool_t AliMagF::LoadParameterization()
 {
   if (fMeasuredMap) {
-    AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
-    return kTRUE;
+    AliFatal(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
   }
   //
   char* fname = gSystem->ExpandPathName(GetDataFileName());
   TFile* file = TFile::Open(fname);
   if (!file) {
-    AliError(Form("Failed to open magnetic field data file %s\n",fname)); 
-    return kFALSE;
+    AliFatal(Form("Failed to open magnetic field data file %s\n",fname)); 
   }
   //
   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
   if (!fMeasuredMap) {
-    AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); 
-    return kFALSE;
+    AliFatal(Form("Did not find field %s in %s\n",GetParamName(),fname)); 
   }
   file->Close();
   delete file;
@@ -239,7 +253,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;
   }
@@ -354,7 +368,7 @@ void AliMagF::MachineField(const Double_t *x, Double_t *b) const
 //_______________________________________________________________________
 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
 {
-  // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
+  // Method to calculate the integral_0^z of br,bt,bz 
   b[0]=b[1]=b[2]=0.0;
   if (fMeasuredMap) {
     fMeasuredMap->GetTPCInt(xyz,b);
@@ -362,10 +376,21 @@ void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
   }
 }
 
+//_______________________________________________________________________
+void AliMagF::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
+{
+  // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
+  b[0]=b[1]=b[2]=0.0;
+  if (fMeasuredMap) {
+    fMeasuredMap->GetTPCRatInt(xyz,b);
+    b[2] /= 100;
+  }
+}
+
 //_______________________________________________________________________
 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 {
-  // Method to calculate the integral of magnetic integral from point to nearest cathode plane
+  // Method to calculate the integral_0^z of br,bt,bz 
   // in cylindrical coordiates ( -pi<phi<pi convention )
   b[0]=b[1]=b[2]=0.0;
   if (fMeasuredMap) {
@@ -374,6 +399,18 @@ void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
   }
 }
 
+//_______________________________________________________________________
+void AliMagF::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
+{
+  // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
+  // in cylindrical coordiates ( -pi<phi<pi convention )
+  b[0]=b[1]=b[2]=0.0;
+  if (fMeasuredMap) {
+    fMeasuredMap->GetTPCRatIntCyl(rphiz,b);
+    b[2] /= 100;
+  }
+}
+
 //_______________________________________________________________________
 void AliMagF::SetFactorSol(Float_t fc)
 {
@@ -437,7 +474,7 @@ AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention,
   const Float_t tolerance=0.03; // relative current tolerance
   const Float_t zero=77.;       // "zero" current (A)
   //
-  BMap_t map;
+  BMap_t map = k5kG;
   double sclL3,sclDip;
   //
   Float_t l3Pol = l3Cur > 0 ? 1:-1;
@@ -449,8 +486,7 @@ AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention,
   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;
+      AliFatalGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
     }
   }
   //
@@ -463,19 +499,17 @@ AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention,
   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)                                { sclL3 = 0;  map  = k5kGUniform;}
+    else if (l3Cur <= zero && diCur<=zero)   { sclL3=0; sclDip=0; map  = k5kGUniform;}
     else {
-      AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
-      return 0;
+      AliFatalGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
     }
   }
   //
-  if (sclDip!=0 && (map==k5kG || map==k2kG) &&
-      ((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 (sclDip!=0 && map!=k5kGUniform) {
+    if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) { 
+      AliFatalGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
+                                    l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
+    }
   }
   //
   if (l3Pol<0) sclL3  = -sclL3;
@@ -496,7 +530,7 @@ AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention,
   // LHC and DCS08 conventions have opposite dipole polarities
   if ( GetPolarityConvention() != convention) sclDip = -sclDip;
   //
-  return new AliMagF("MagneticFieldMap", ttl, 2, sclL3, sclDip, 10., map, path,btype,beamenergy);
+  return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
   //
 }
 
@@ -514,3 +548,19 @@ const char*  AliMagF::GetBeamTypeText() const
   }
 }
 
+//_____________________________________________________________________________
+void AliMagF::Print(Option_t *opt) const
+{
+  // print short or long info
+  TString opts = opt; opts.ToLower();
+  AliInfo(Form("%s:%s",GetName(),GetTitle()));
+  AliInfo(Form("Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
+              GetFactorSol(),(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
+              fDipoleOFF ? "OFF":"ON",GetFactorDip(),fMapType==k5kGUniform?" |Constant Field!":""));
+  if (opts.Contains("a")) {
+    AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
+                fBeamType==kBeamTypeAA ? "A-A":(fBeamType==kBeamTypepp ? "p-p":"OFF"),
+                fBeamEnergy,fQuadGradient,fDipoleField));
+    AliInfo(Form("Uses %s of %s",GetParamName(),GetDataFileName()));
+  }
+}