Magnetic field map for ALICE for L3+muon spectrometer stored in 3 seperate
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Feb 2002 11:41:28 +0000 (11:41 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Feb 2002 11:41:28 +0000 (11:41 +0000)
root files.

STEER/AliFieldMap.cxx [new file with mode: 0644]
STEER/AliFieldMap.h [new file with mode: 0644]
STEER/AliMagFMaps.cxx [new file with mode: 0644]
STEER/AliMagFMaps.h [new file with mode: 0644]

diff --git a/STEER/AliFieldMap.cxx b/STEER/AliFieldMap.cxx
new file mode 100644 (file)
index 0000000..0e54c4a
--- /dev/null
@@ -0,0 +1,222 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+//
+// Author: Andreas Morsch <andreas.morsch@cern.ch>
+//
+
+#include <TVector.h>
+#include "AliFieldMap.h"
+#include "TSystem.h"
+
+ClassImp(AliFieldMap)
+
+//________________________________________
+AliFieldMap::AliFieldMap()
+{
+  //
+  // Standard constructor
+  //
+  fB = 0;
+}
+
+AliFieldMap::AliFieldMap(const char *name, const char *title)
+    : TNamed(name,title)
+{
+  //
+  // Standard constructor
+  //
+  fB = 0;
+  ReadField();
+}
+
+AliFieldMap::~AliFieldMap()
+{
+//
+// Destructor
+//  
+  delete fB;
+}
+
+//________________________________________
+AliFieldMap::AliFieldMap(const AliFieldMap &map)
+{
+  //
+  // Copy constructor
+  //
+  map.Copy(*this);
+}
+
+//________________________________________
+void AliFieldMap::ReadField()
+{
+  // 
+  // Method to read the magnetic field map from file
+  //
+  FILE* magfile;
+  FILE* endf = fopen("end.table", "r");
+  FILE* out  = fopen("out", "w");
+  
+  Int_t   ix, iy, iz, ipx, ipy, ipz;
+  Float_t bx, by, bz;
+  char *fname = 0;
+  printf("%s: Reading Magnetic Field Map %s from file %s\n",
+        ClassName(),fName.Data(),fTitle.Data());
+
+  fname   = gSystem->ExpandPathName(fTitle.Data());
+  magfile = fopen(fname,"r");
+  delete [] fname;
+
+  fscanf(magfile,"%d %d %d %f %f %f %f %f %f", 
+        &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
+  fXdeli = 1./fXdel;
+  fYdeli = 1./fYdel;   
+  fZdeli = 1./fZdel;   
+  fXend  = fXbeg + (fXn-1)*fXdel;
+  fYend  = fYbeg + (fYn-1)*fYdel;
+  fZend  = fZbeg + (fZn-1)*fZdel;
+  
+  Int_t nDim   = fXn*fYn*fZn;
+
+//  Float_t x,y,z,b;
+
+  fB = new TVector(3*nDim);
+  if (magfile) {
+      for (ix = 0; ix < fXn; ix++) {
+         ipx=ix*3*(fZn*fYn);
+         for (iy = 0; iy < fYn; iy++) {
+             ipy=ipx+iy*3*fZn;
+             for (iz = 0; iz < fZn; iz++) {
+                 ipz=ipy+iz*3;
+
+                 if (iz == -1) 
+                     fscanf(endf,"%f %f %f", &bx,&by,&bz);
+                 else if (iz > -1)
+                     fscanf(magfile," %f %f %f", &bx, &by, &bz);
+                 else 
+                     continue;
+
+//               fscanf(magfile,"%f %f %f %f %f %f %f ",
+//                      &x, &y, &z, &bx,&by,&bz, &b);
+//               fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
+
+                 (*fB)(ipz+2) = 10.*bz;
+                 (*fB)(ipz+1) = 10.*by;
+                 (*fB)(ipz  ) = 10.*bx;
+             } //iz
+         } // iy
+      } // ix
+/*
+//
+// this part for interpolation between z = 700 and 720 cm to get
+// z = 710 cm
+//      
+      for (ix = 0; ix < fXn; ix++) {
+         ipx=ix*3*(fZn*fYn);
+         for (iy = 0; iy < fYn; iy++) {
+             ipy=ipx+iy*3*fZn;
+             Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
+             Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
+             Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;           
+             ipz=ipy+3;
+             (*fB)(ipz+2) = bzz;
+             (*fB)(ipz+1) = byy;
+             (*fB)(ipz  ) = bxx;
+         } // iy
+      } // ix
+*/      
+  } else { 
+      printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
+      exit(1);
+  } // if mafile
+}
+
+void AliFieldMap::Field(Float_t *x, Float_t *b)
+{
+//
+// Use simple interpolation to obtain field at point x
+//
+    Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
+       bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
+    const Double_t kone=1;
+    Int_t ix, iy, iz;
+    b[0]=b[1]=b[2]=0;
+//
+    
+    xl[0]=TMath::Abs(x[0])-fXbeg;
+    xl[1]=TMath::Abs(x[1])-fYbeg;
+    xl[2]=x[2]-fZbeg;
+    
+    hix=xl[0]*fXdeli;
+    ratx=hix-int(hix);
+    ix=int(hix);
+    
+    hiy=xl[1]*fYdeli;
+    raty=hiy-int(hiy);
+    iy=int(hiy);
+    
+    hiz=xl[2]*fZdeli;
+    ratz=hiz-int(hiz);
+    iz=int(hiz);
+
+    ratx1=kone-ratx;
+    raty1=kone-raty;
+    ratz1=kone-ratz;
+
+    bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
+    bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
+    blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
+    blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
+    bhz   = blyhz             *raty1+bhyhz             *raty;
+    blz   = blylz             *raty1+bhylz             *raty;
+    b[0]  = blz               *ratz1+bhz               *ratz;
+    //
+    bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
+    bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
+    blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
+    blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
+    bhz   = blyhz             *raty1+bhyhz             *raty;
+    blz   = blylz             *raty1+bhylz             *raty;
+    b[1]  = blz               *ratz1+bhz               *ratz;
+    //
+    bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
+    bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
+    blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
+    blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
+    bhz   = blyhz             *raty1+bhyhz             *raty;
+    blz   = blylz             *raty1+bhylz             *raty;
+    b[2]  = blz               *ratz1+bhz               *ratz;
+}
+
+//________________________________________
+void AliFieldMap::Copy(AliFieldMap & /* magf */) const
+{
+  //
+  // Copy *this onto magf -- Not implemented
+  //
+  Fatal("Copy","Not implemented!\n");
+}
+
+//________________________________________
+AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
+{
+  magf.Copy(*this);
+  return *this;
+}
diff --git a/STEER/AliFieldMap.h b/STEER/AliFieldMap.h
new file mode 100644 (file)
index 0000000..0f6cac2
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALIFIELDMAP_H
+#define ALIFIELDMAP_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//
+// Author: Andreas Morsch <andreas.morsch@cern.ch>
+//
+#include <TNamed.h>
+#include <TVector.h>
+
+class AliFieldMap : public TNamed
+{
+  //Alice Magnetic Field with constant mesh
+
+public:
+    AliFieldMap();
+    AliFieldMap(const char *name, const char *title);
+    AliFieldMap(const AliFieldMap &mag);
+    virtual ~AliFieldMap();
+    void Copy(AliFieldMap &map) const;
+    virtual AliFieldMap & operator=(const AliFieldMap &map);
+
+    virtual void Field(Float_t *x, Float_t *b);
+    Float_t Bx(const Int_t ix, const Int_t iy, const Int_t iz) {
+       return (*fB)(3*(ix*(fZn*fYn)+iy*fZn+iz));
+    }
+    Float_t By(const Int_t ix, const Int_t iy, const Int_t iz) {
+       return (*fB)(3*(ix*(fZn*fYn)+iy*fZn+iz)+1);
+    }
+    Float_t Bz(const Int_t ix, const Int_t iy, const Int_t iz) {
+       return (*fB)(3*(ix*(fZn*fYn)+iy*fZn+iz)+2);
+    }
+
+    Bool_t Inside(Float_t x, Float_t y, Float_t z) 
+       { return (x > fXbeg && x <= fXend &&
+                 y > fYbeg && y <= fYend &&
+                 z > fZbeg && z <= fZend);
+       }
+    Float_t Xmin()  {return fXbeg;}
+    Float_t Xmax()  {return fXend;}
+    Float_t DelX()  {return fXdel;}
+    Float_t DeliX() {return fXdeli;}
+    
+    Float_t Ymin()  {return fYbeg;}
+    Float_t Ymax()  {return fYend;}
+    Float_t DelY()  {return fYdel;}
+    Float_t DeliY() {return fYdeli;}
+    
+    Float_t Zmin()  {return fZbeg;}
+    Float_t Zmax()  {return fZend;}
+    Float_t DelZ()  {return fZdel;}
+    Float_t DeliZ() {return fZdeli;}
+ private:
+    void    ReadField();
+ protected:
+
+    Float_t    fXbeg;     // Start of mesh in x
+    Float_t    fYbeg;     // Start of mesh in y
+    Float_t    fZbeg;     // Start of mesh in z
+    Float_t    fXend;     // End of mesh in x
+    Float_t    fYend;     // End of mesh in y
+    Float_t    fZend;     // End of mesh in z
+    Float_t    fXdel;     // Mesh step in x
+    Float_t    fYdel;     // Mesh step in y
+    Float_t    fZdel;     // Mesh step in z
+    Double_t   fXdeli;    // Inverse of Mesh step in x
+    Double_t   fYdeli;    // Inverse of Mesh step in y
+    Double_t   fZdeli;    // Inverse of Mesh step in z
+    Int_t      fXn;       // Number of mesh points in x
+    Int_t      fYn;       // Number of mesh points in y
+    Int_t      fZn;       // Number of mesh points in z
+    TVector*   fB;        // Field map
+    
+    ClassDef(AliFieldMap,1)  //Class for Field Map
+};
+
+#endif
diff --git a/STEER/AliMagFMaps.cxx b/STEER/AliMagFMaps.cxx
new file mode 100644 (file)
index 0000000..e2ecef0
--- /dev/null
@@ -0,0 +1,243 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+//
+// Author: Andreas Morsch <andreas.morsch@cern.ch>
+//
+
+#include <TFile.h>
+#include <TSystem.h>
+#include "AliFieldMap.h"
+#include "AliMagFMaps.h"
+
+
+ClassImp(AliMagFMaps)
+
+//________________________________________
+AliMagFMaps::AliMagFMaps(const char *name, const char *title, const Int_t integ, 
+                    const Float_t factor, const Float_t fmax, const Int_t map)
+  : AliMagF(name,title,integ,factor,fmax)
+{
+  //
+  // Standard constructor
+  //
+  fType         = kConMesh;
+  fFieldMap[0]  = 0;
+  char* fname;
+  
+  fMap = map;
+  TFile* file = 0;
+  if(fDebug>-1) printf("%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
+        ClassName(),fName.Data(), fMap, factor,fTitle.Data());
+    if (fMap == k2kG) {
+/*
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/L3B02.root");
+       file = new TFile(fname);
+       fFieldMap[0] = (AliFieldMap*) file->Get("L3B02");
+       file->Close();
+       delete file;
+
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/DipB02.root");
+       file = new TFile(fname);
+       fFieldMap[1] = (AliFieldMap*) file->Get("DipB02");
+       file->Close();
+       delete file;;
+
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/ExtB02.root");
+       file = new TFile(fname);
+       fFieldMap[2] = (AliFieldMap*) file->Get("ExtB02");
+       file->Close();
+       delete file;
+*/
+       Fatal("AliMagFMaps", "Bcent = 0.2 T not yet available");
+       
+       fSolenoid = 2.;
+    } else if (fMap == k4kG) {
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/L3B04.root");
+       file = new TFile(fname);
+       fFieldMap[0] = (AliFieldMap*) file->Get("L3B04");
+       file->Close();
+       delete file;
+
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/DipB04.root");
+       file = new TFile(fname);
+       fFieldMap[1] = (AliFieldMap*) file->Get("DipB04");
+       file->Close();
+       delete file;;
+
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/ExtB04.root");
+       file = new TFile(fname);
+       fFieldMap[2] = (AliFieldMap*) file->Get("ExtB04");
+       file->Close();
+       delete file;
+       fSolenoid = 4.;
+    } else if (fMap == k5kG) {
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/L3B05.root");
+       file = new TFile(fname);
+       fFieldMap[0] = (AliFieldMap*) file->Get("L3B05");
+       file->Close();
+       delete file;
+
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/DipB05.root");
+       file = new TFile(fname);
+       fFieldMap[1] = (AliFieldMap*) file->Get("DipB05");
+       file->Close();
+       delete file;;
+
+       fname = gSystem->ExpandPathName("$(ALICE_ROOT)/data/maps/ExtB05.root");
+       file = new TFile(fname);
+       fFieldMap[2] = (AliFieldMap*) file->Get("ExtB05");
+       file->Close();
+       delete file;
+
+       fSolenoid = 5.;
+    }
+    SetL3ConstField(0);
+}
+
+//________________________________________
+AliMagFMaps::AliMagFMaps(const AliMagFMaps &magf)
+{
+  //
+  // Copy constructor
+  //
+  magf.Copy(*this);
+}
+
+AliMagFMaps::~AliMagFMaps()
+{
+//
+//  Destructor
+//
+    delete fFieldMap[0];
+    delete fFieldMap[1];
+    delete fFieldMap[2];    
+}
+
+
+Float_t AliMagFMaps::SolenoidField() const
+{
+//
+// Returns max. L3 (solenoid) field strength 
+// according to field map setting
+
+    return fSolenoid;
+}
+
+    
+
+//________________________________________
+void AliMagFMaps::Field(Float_t *x, Float_t *b)
+{
+  //
+  // Method to calculate the magnetic field
+  //
+  Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
+    bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
+  const Double_t kone=1;
+  Int_t ix, iy, iz;
+  
+  // --- find the position in the grid ---
+  
+  b[0]=b[1]=b[2]=0;
+  AliFieldMap* map = 0;
+
+  if (fFieldMap[0]->Inside(x[0], x[1], x[2])) {
+      map = fFieldMap[0];
+      if (fL3Option) {
+//
+//     Constant L3 field, if this option was selected
+//
+         b[2] = fSolenoid;
+         return;
+      }
+  } else if (fFieldMap[1]->Inside(x[0], x[1], x[2])) {
+      map = fFieldMap[1];
+  } else if (fFieldMap[2]->Inside(x[0], x[1], x[2])) {
+      map = fFieldMap[2];
+  }
+  
+  if(map){
+      map->Field(x,b);
+  } else {
+//This is the ZDC part
+      Float_t rad2=x[0]*x[0]+x[1]*x[1];
+      if(x[2]>kCORBEG2 && x[2]<kCOREND2){
+         if(rad2<kCOR2RA2){
+             b[0] = kFCORN2;
+         }
+      }
+      else if(x[2]>kZ1BEG && x[2]<kZ1END){  
+         if(rad2<kZ1RA2){
+             b[0] = -kG1*x[1];
+             b[1] = -kG1*x[0];
+         }
+      }
+      else if(x[2]>kZ2BEG && x[2]<kZ2END){  
+         if(rad2<kZ2RA2){
+             b[0] = kG1*x[1];
+             b[1] = kG1*x[0];
+         }
+      }
+      else if(x[2]>kZ3BEG && x[2]<kZ3END){  
+         if(rad2<kZ3RA2){
+             b[0] = kG1*x[1];
+             b[1] = kG1*x[0];
+         }
+      }
+      else if(x[2]>kZ4BEG && x[2]<kZ4END){  
+         if(rad2<kZ4RA2){
+             b[0] = -kG1*x[1];
+             b[1] = -kG1*x[0];
+         }
+      }
+      else if(x[2]>kD1BEG && x[2]<kD1END){ 
+         if(rad2<kD1RA2){
+             b[1] = -kFDIP;
+         }
+      }
+      else if(x[2]>kD2BEG && x[2]<kD2END){
+         if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
+            || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
+             b[1] = kFDIP;
+         }
+      }
+  }
+  if(fFactor!=1) {
+      b[0]*=fFactor;
+      b[1]*=fFactor;
+      b[2]*=fFactor;
+  }
+}
+
+//________________________________________
+void AliMagFMaps::Copy(AliMagFMaps & /* magf */) const
+{
+  //
+  // Copy *this onto magf -- Not implemented
+  //
+  Fatal("Copy","Not implemented!\n");
+}
+
+//________________________________________
+AliMagFMaps & AliMagFMaps::operator =(const AliMagFMaps &magf)
+{
+  magf.Copy(*this);
+  return *this;
+}
diff --git a/STEER/AliMagFMaps.h b/STEER/AliMagFMaps.h
new file mode 100644 (file)
index 0000000..262ac7d
--- /dev/null
@@ -0,0 +1,40 @@
+#ifndef ALIMAGFMAPS_H
+#define ALIMAGFMAPS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//
+// Author: Andreas Morsch <andreas.morsch@cern.ch>
+//
+
+#include "AliMagF.h"
+class AliFieldMap;
+
+class AliMagFMaps : public AliMagF
+{
+  //Alice Magnetic Field with constant mesh
+
+public:
+    enum constants {k2kG, k4kG, k5kG};
+    AliMagFMaps(){fFieldMap[0] = fFieldMap[1] = fFieldMap[2] = 0;}
+    AliMagFMaps(const char *name, const char *title, const Int_t integ,
+               const Float_t factor, const Float_t fmax, const Int_t map = k2kG);
+    AliMagFMaps(const AliMagFMaps &mag);
+    virtual ~AliMagFMaps();
+    virtual void    Field(Float_t *x, Float_t *b);
+    AliFieldMap* FieldMap(Int_t i) {return fFieldMap[i];}
+    virtual Float_t SolenoidField() const;
+    virtual void    SetL3ConstField(Int_t flag = 0) {fL3Option = flag;}
+    
+    void Copy(AliMagFMaps &magf) const;
+    virtual AliMagFMaps & operator=(const AliMagFMaps &magf);
+protected:
+  AliFieldMap* fFieldMap[3];     // Field maps
+  Float_t      fSolenoid;        // Solenoid field setting
+  Int_t        fL3Option;        // Option for field inside L3
+  ClassDef(AliMagFMaps,1)        // Class for all Alice MagField using three Maps with Constant Mesh
+};
+
+#endif