PWGLF/FORWARD/analysis2/sim/GRP.C

   1 /**
   2  * @file   GRP.C
   3  * @author Christian Holm Christensen <cholm@nbi.dk>
   4  * @date   Wed Aug 20 10:00:39 2014
   5  *
   6  * @brief  Class that hold summary of GRP data.
   7  *
   8  *
   9  */
  10 // #include <fstream>
  11
  12 /**
  13  * Class that hold summary of GRP data, and has methods to obtain
  14  * these, either from a previously dumped file, or directly from OCDB.
  15  *
  16  * This is used by the simulation setup to ensure we have the proper
  17  * beam settings for an anchor run.
  18  */
  19 struct GRPData
  20 {
  21   /**
  22    * A beam
  23    */
  24   struct Beam
  25   {
  26     UInt_t z;
  27     UInt_t a;
  28     /**
  29      * Get the per nucleon beam energy given a per charge beam energy
  30      *
  31      * @param e Per charge beam energy
  32      *
  33      * @return Per nucleon beam energy
  34      */
  35     Float_t GetEnergy(Float_t e) const
  36     {
  37       return e * sqrt(float(z)/a);
  38     }
  39     /**
  40      * Check if this is a proton
  41      *
  42      * @return True if z and a are 1
  43      */
  44     Bool_t IsP() const { return z == 1 && a == 1; }
  45     /**
  46      * Check if this is a nucleus
  47      *
  48      * @return true if not P
  49      */
  50     Bool_t IsA() const { return !IsP(); }
  51     /**
  52      * Get the beam type name
  53      *
  54      * @return P or A
  55      */
  56     const char* Name() { return (IsP() ? "P" : IsA() ? "A" : ""); }
  57     /**
  58      * Set the beam type from single beam spec from LHC
  59      *
  60      * @param b Beam type from LHC
  61      */
  62     void Set(UInt_t b)
  63     {
  64       z = b % 1000;
  65       a = b / 1000;
  66     }
  67     /**
  68      * Set from either a full LHC spec or from ALICE spec
  69      *
  70      * @param b Beam
  71      */
  72     void Set(const TString& b)
  73     {
  74       b.ToUpper();
  75       a = 0;
  76       z = 0;
  77       if (b.EqualTo("P") || b.EqualTo("PROTON")) {
  78         z = a = 1;
  79       }
  80       else if (b.EqualTo("A") || b.BeginsWidth("PB")) {
  81         z = 82;
  82         a = 208;
  83       }
  84     }
  85   };
  86   UInt_t  beamEnergy; // Total energy in center of mass
  87   UInt_t  energy; // Center of mass energy per nucleon
  88   TString period; // The period
  89   UInt_t  run;   // The run number
  90   Beam    beam1; // Target beam
  91   Beam    beam2; // Projectile beam
  92   /**
  93    * Constructor.
  94    *
  95    * @param r Run number
  96    */
  97   GRPData(UInt_t r)
  98     : beamEnergy(0), energy(0), period(""), beam1(), beam2(), run(r)
  99   {
 100     Info("GRP", "Will try from file %s", FileName());
 101     if (!ReadFromFile()) {
 102       Warning("GRP", "Failed to read from file, trying from OCDB");
 103       if (!ReadFromOCDB(run))
 104         Fatal("GRP", "Failed to get GRP values");
 105     }
 106     Print();
 107   }
 108   /**
 109    * Print information
 110    *
 111    */
 112   void Print()
 113   {
 114     Printf("%s run %d\n"
 115            "  Beam energy: %d GeV\n"
 116            "  sqrt(s_NN):  %d GeV\n"
 117            "  Beam 1:      %s (%d %d)\n"
 118            "  Beam 2:      %s (%d %d)",
 119            period.Data(), run, beamEnergy, energy,
 120            beam1.Name(), beam1.a, beam1.z,
 121            beam2.Name(), beam2.a, beam2.z);
 122   }
 123   /**
 124    * Check if this run was a p-p run
 125    *
 126    * @return Return true if both beams are p
 127    */
 128   Bool_t IsPP() const { return beam1.IsP() && beam2.IsP(); }
 129   /**
 130    * Check if this run was a p-A run
 131    *
 132    * @return Return true beam 1 is p and 2 is A
 133    */
 134   Bool_t IsPA() const { return beam1.IsP() && beam2.IsA(); }
 135   /**
 136    * Check if this run was a A-p run
 137    *
 138    * @return Return true beam 1 is A and 2 is p
 139    */
 140   Bool_t IsAP() const { return beam1.IsA() && beam2.IsP(); }
 141   /**
 142    * Check if this run was a A-A run
 143    *
 144    * @return Return true if both beams are A
 145    */
 146   Bool_t IsAA() const { return beam1.IsA() && beam2.IsA(); }
 147   /**
 148    * Get the year
 149    *
 150    * @return Year
 151    */
 152   UInt_t Year() const
 153   {
 154     TString tmp = period(3,2);
 155     return tmp.Atoi() + 2000;
 156   }
 157   const char* FileName() const { return "grp.dat"; }
 158   /**
 159    * Get unsigned int value from string value of TObjString
 160    *
 161    * @param o Object
 162    *
 163    * @return value
 164    */
 165   static UInt_t Obj2UInt(const TObject* o)
 166   {
 167     return (static_cast<TObjString*>(o))->String().Atoi();
 168   }
 169   /**
 170    * Read values from a file
 171    *
 172    * @return true on success
 173    */
 174   Bool_t ReadFromFile()
 175   {
 176     if (gSystem->AccessPathName(FileName())) {
 177       Info("GRP", "Cannot open file %s", FileName());
 178       return false;
 179     }
 180
 181     std::ifstream* pin = new std::ifstream("grp.dat");
 182     std::ifstream& in  = *pin;
 183     TString line;
 184     TString env;
 185     do {
 186       line.ReadLine(in);
 187       if (line.IsNull()) continue;
 188       if (line[0] == '#') continue;
 189       env = line;
 190       break;
 191     } while (!in.eof());
 192     if (env.IsNull()) {
 193       Info("GRP", "No line found");
 194       return false;
 195     }
 196     Info("GRP", "Got the line %s", env.Data());
 197     TObjArray* tokens = env.Tokenize(":");
 198     if (tokens->GetEntries() < 8) {
 199       Warning("GRP", "Failed to get enough field from GRP_SUMMARY=%s",
 200               env.Data());
 201       return false;
 202     }
 203     period     = tokens->At(0)->GetName();
 204     run        = Obj2UInt(tokens->At(1));
 205     beamEnergy = Obj2UInt(tokens->At(2));
 206     energy     = Obj2UInt(tokens->At(3));
 207     beam1.a    = Obj2UInt(tokens->At(4));
 208     beam1.z    = Obj2UInt(tokens->At(5));
 209     beam2.a    = Obj2UInt(tokens->At(6));
 210     beam2.z    = Obj2UInt(tokens->At(7));
 211     pin->close();
 212     delete pin;
 213     return true;
 214   }
 215   /**
 216    * Read values from OCDB.  Also dumps values to file.
 217    *
 218    * @param r run number
 219    *
 220    * @return true on success
 221    */
 222   Bool_t ReadFromOCDB(UInt_t r)
 223   {
 224      AliCDBManager* cdb = AliCDBManager::Instance();
 225      cdb->SetDefaultStorageFromRun(r);
 226      cdb->SetRun(r);
 227      AliCDBEntry*   ent = cdb->Get("GRP/GRP/Data");
 228      if (!ent) {
 229        Warning("GRP", "No GRP data found");
 230        cdb->SetRun(-1);
 231        return false;
 232      }
 233      // Info("GRP", "Got GRP");
 234      // ent->PrintMetaData();
 235
 236      AliGRPObject*  obj        = static_cast<AliGRPObject*>(ent->GetObject());
 237      obj->Print();
 238      run                       = r;
 239      period                    = obj->GetLHCPeriod();
 240      beamEnergy                = TMath::Ceil(obj->GetBeamEnergy());
 241      TString        beam1T     = obj->GetSingleBeamType(0);
 242      TString        beam2T     = obj->GetSingleBeamType(1);
 243
 244      if (!beam1T.IsNull() && !beam2T.IsNull()) {
 245        beam1.Set(beam1T.Atoi());
 246        beam2.Set(beam2T.Atoi());
 247      }
 248      else {
 249        TString beamType = obj->GetBeamType();
 250        if (beamType.IsNull()) {
 251          Warning("GRP", "No beam type");
 252          cdb->SetRun(-1);
 253          return false;
 254        }
 255
 256        TObjArray* ab  = beamType.Tokenize("-");
 257        ab->ls();
 258        beam1T = ab->At(0)->GetName();
 259        beam2T = ab->At(1)->GetName();
 260        beam1.Set(beam1T);
 261        beam2.Set(beam2T);
 262      }
 263      // Massage the beam energy in case we had sqrt{s_NN} instead of
 264      // beam energy.
 265      if (TMath::Abs(beamEnergy - 1379.8) < 10 && beam1.IsA() && beam2.IsA())
 266        beamEnergy = 3500;
 267      energy = TMath::Ceil(2*beamEnergy * TMath::Sqrt(Float_t(beam1.z*beam2.z)/
 268                                                      (beam1.a*beam2.a)));
 269
 270      const char* fn = FileName();
 271      std::ofstream* pout = new std::ofstream(fn);
 272      std::ofstream& out  = *pout;
 273      out << "# PERIOD:RUN:BEAMENERGY:ENERGY:BEAM1A:BEAM1Z:BEAM2A:BEAM2Z\n"
 274          << period     << ":"
 275          << run        << ":"
 276          << beamEnergy << ":"
 277          << energy     << ":"
 278          << beam1.a    << ":"
 279          << beam1.z    << ":"
 280          << beam2.a    << ":"
 281          << beam2.z    << "\n"
 282          << "# EOF "   << std::endl;
 283      pout->close();
 284      delete pout;
 285      cdb->SetRun(-1);
 286      return true;
 287   }
 288   /**
 289    * Get the value of an environment variable as a unsigned int
 290    *
 291    * @param envName Enviroment variable name
 292    * @param def     Default value
 293    *
 294    * @return As unsigned int, or default
 295    */
 296   static UInt_t Env2UInt(const char* envName, UInt_t def)
 297   {
 298     TString val(gSystem->Getenv(envName));
 299     if (val.IsNull()) return def;
 300     return UInt_t(val.Atoll());
 301   }
 302 };
 303 GRPData* grp;
 304 void GRP(UInt_t run=0)
 305 {
 306   grp = new GRPData(run);
 307 }
 308 //
 309 // EOF
 310 //
 311
 312
 313
 314