]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliSimulation.cxx
Set beam energy in GRP object to energy per charge.
[u/mrichter/AliRoot.git] / STEER / AliSimulation.cxx
index b6fff7bcde399a1be34b92995d4f5dce8cdfd936..ca5b3a6948c6a5f793618bbb2a89ef1ca7009a8b 100644 (file)
 #include <TSystem.h>
 #include <TVirtualMC.h>
 #include <TVirtualMCApplication.h>
+#include <TDatime.h>
 
 #include "AliAlignObj.h"
 #include "AliCDBEntry.h"
 #include "AliCDBManager.h"
+#include "AliGRPManager.h"
 #include "AliCDBStorage.h"
 #include "AliCTPRawData.h"
 #include "AliCentralTrigger.h"
@@ -165,7 +167,7 @@ AliSimulation::AliSimulation(const char* configFileName,
   fLoadAlObjsListOfDets("ALL"),
   fMakeSDigits("ALL"),
   fMakeDigits("ALL"),
-  fMakeTrigger(""),
+  fTriggerConfig(""),
   fMakeDigitsFromHits(""),
   fWriteRawData(""),
   fRawDataFileName(""),
@@ -190,11 +192,16 @@ AliSimulation::AliSimulation(const char* configFileName,
   fSetRunNumberFromDataCalled(kFALSE),
   fEmbeddingFlag(kFALSE),
   fLego(NULL),
+  fKey(0),
+  fUseVertexFromCDB(0),
+  fUseMagFieldFromGRP(0),
+  fGRPWriteLocation(Form("local://%s", gSystem->pwd())),
   fQADetectors("ALL"),                  
   fQATasks("ALL"),     
   fRunQA(kTRUE), 
   fEventSpecie(AliRecoParam::kDefault),
   fWriteQAExpertData(kTRUE), 
+  fGeometryFile(),
   fRunHLT("default"),
   fpHLT(NULL),
   fWriteGRPEntry(kTRUE)
@@ -360,7 +367,8 @@ void AliSimulation::SetCDBLock() {
   // Set CDB lock: from now on it is forbidden to reset the run number
   // or the default storage or to activate any further storage!
   
-  AliCDBManager::Instance()->SetLock(1);
+  ULong_t key = AliCDBManager::Instance()->SetLock(1);
+  if (key) fKey = key;
 }
 
 //_____________________________________________________________________________
@@ -525,7 +533,7 @@ Bool_t AliSimulation::MisalignGeometry(AliRunLoader *runLoader)
   }
   
   // Export ideal geometry 
-  if(!gAlice->IsRootGeometry()) AliGeomManager::GetGeometry()->Export("geometry.root");
+  if(!IsGeometryFromFile()) AliGeomManager::GetGeometry()->Export("geometry.root");
 
   // Load alignment data from CDB and apply to geometry through AliGeomManager
   if(fLoadAlignFromCDB){
@@ -601,7 +609,7 @@ Bool_t AliSimulation::Run(Int_t nEvents)
 // run the generation, simulation and digitization
 
  
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
   AliSysInfo::AddStamp("Start_Run");
   
   // Load run number and seed from environmental vars
@@ -643,9 +651,7 @@ Bool_t AliSimulation::Run(Int_t nEvents)
   if (!AliGeomManager::GetGeometry()) {
     // Initialize the geometry manager
     AliGeomManager::LoadGeometry("geometry.root");
-  AliSysInfo::AddStamp("GetGeometry");
-
-    
+    AliSysInfo::AddStamp("GetGeometry");
 //    // Check that the consistency of symbolic names for the activated subdetectors
 //    // in the geometry loaded by AliGeomManager
 //    AliRunLoader* runLoader = LoadRun("READ");
@@ -706,7 +712,7 @@ Bool_t AliSimulation::Run(Int_t nEvents)
   
   
   // digits -> trigger
-  if (!RunTrigger(fMakeTrigger,fMakeDigits)) {
+  if (!RunTrigger(fTriggerConfig,fMakeDigits)) {
     if (fStopOnError) return kFALSE;
   }
 
@@ -794,7 +800,7 @@ Bool_t AliSimulation::RunLego(const char *setup, Int_t nc1, Float_t c1min,
 
 // run the generation and simulation
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
 
   // initialize CDB storage and run number from external environment
   // (either CDB manager or AliSimulation setters)
@@ -851,10 +857,10 @@ Bool_t AliSimulation::RunLego(const char *setup, Int_t nc1, Float_t c1min,
   
   AliPDG::AddParticlesToPdgDataBase();  
   
-  gAlice->GetMCApp()->Init();
-  
   gMC->SetMagField(TGeoGlobalMagField::Instance()->GetField());
 
+  gAlice->GetMCApp()->Init();
+  
   
   //Must be here because some MCs (G4) adds detectors here and not in Config.C
   gAlice->InitLoaders();
@@ -892,7 +898,7 @@ Bool_t AliSimulation::RunTrigger(const char* config, const char* detectors)
 {
   // run the trigger
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
 
   // initialize CDB storage from external environment
   // (either CDB manager or AliSimulation setters),
@@ -912,8 +918,8 @@ Bool_t AliSimulation::RunTrigger(const char* config, const char* detectors)
    TString trconfiguration = config;
 
    if (trconfiguration.IsNull()) {
-     if (strcmp(gAlice->GetTriggerDescriptor(),"")) {
-       trconfiguration = gAlice->GetTriggerDescriptor();
+     if(!fTriggerConfig.IsNull()) {
+       trconfiguration = fTriggerConfig;
      }
      else
        AliWarning("No trigger descriptor is specified. Loading the one that is in the CDB.");
@@ -955,7 +961,7 @@ Bool_t AliSimulation::RunSimulation(Int_t nEvents)
 {
 // run the generation and simulation
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
 
   // initialize CDB storage and run number from external environment
   // (either CDB manager or AliSimulation setters)
@@ -981,13 +987,50 @@ Bool_t AliSimulation::RunSimulation(Int_t nEvents)
   //
   // Initialize ALICE Simulation run
   //
-
   gAlice->Announce();
 
+  //
+  // If requested set the mag. field from the GRP entry.
+  // After this the field is loccked and cannot be changed by Config.C
+  if (fUseMagFieldFromGRP) {
+      AliGRPManager grpM;
+      grpM.ReadGRPEntry();
+      grpM.SetMagField();
+      AliInfo("Field is locked now. It cannot be changed in Config.C");
+  }
+//
+// Execute Config.C
   gROOT->LoadMacro(fConfigFileName.Data());
   gInterpreter->ProcessLine(gAlice->GetConfigFunction());
   AliSysInfo::AddStamp("RunSimulation_Config");
 
+//
+// If requested obtain the vertex position and vertex sigma_z from the CDB
+// This overwrites the settings from the Config.C  
+  if (fUseVertexFromCDB) {
+      Double_t vtxPos[3] = {0., 0., 0.}; 
+      Double_t vtxSig[3] = {0., 0., 0.};
+      AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/Calib/MeanVertex");
+      AliESDVertex* vertex = dynamic_cast<AliESDVertex*> (entry->GetObject());
+      if(vertex) {
+         if(vertex->GetXRes()>2.8) { // > pipe radius --> it's a dummy object, don't use it 
+             entry = AliCDBManager::Instance()->Get("GRP/Calib/MeanVertexSPD");
+             vertex = dynamic_cast<AliESDVertex*> (entry->GetObject());
+         }
+      }
+      if (vertex) {
+         vertex->GetXYZ(vtxPos);
+         vertex->GetSigmaXYZ(vtxSig);
+         AliInfo("Overwriting Config.C vertex settings !");
+         AliInfo(Form("Vertex position from OCDB entry: x = %13.3f, y = %13.3f, z = %13.3f (sigma = %13.3f)\n",
+                vtxPos[0], vtxPos[1], vtxPos[2], vtxSig[2]));
+
+         AliGenerator *gen = gAlice->GetMCApp()->Generator();
+         gen->SetOrigin(vtxPos[0], vtxPos[1], vtxPos[2]);   // vertex position
+         gen->SetSigmaZ(vtxSig[2]);
+      }
+  }
+  
   if(AliCDBManager::Instance()->GetRun() >= 0) { 
     AliRunLoader::Instance()->SetRunNumber(AliCDBManager::Instance()->GetRun());
     AliRunLoader::Instance()->SetNumberOfEventsPerRun(fNEvents);
@@ -996,15 +1039,16 @@ Bool_t AliSimulation::RunSimulation(Int_t nEvents)
   }
   
    AliRunLoader::Instance()->CdGAFile();
-    
-   AliPDG::AddParticlesToPdgDataBase();  
+   
 
-   gAlice->GetMCApp()->Init();
-   AliSysInfo::AddStamp("RunSimulation_InitMCApp");
+   AliPDG::AddParticlesToPdgDataBase();  
 
    gMC->SetMagField(TGeoGlobalMagField::Instance()->GetField());
    AliSysInfo::AddStamp("RunSimulation_GetField");
    
+   gAlice->GetMCApp()->Init();
+   AliSysInfo::AddStamp("RunSimulation_InitMCApp");
+
    //Must be here because some MCs (G4) adds detectors here and not in Config.C
    gAlice->InitLoaders();
    AliRunLoader::Instance()->MakeTree("E");
@@ -1019,15 +1063,6 @@ Bool_t AliSimulation::RunSimulation(Int_t nEvents)
    AliSysInfo::AddStamp("RunSimulation_InitLoaders");
   //___________________________________________________________________________________________
   
-  // Get the trigger descriptor string
-  // Either from AliSimulation or from
-  // gAlice
-  if (fMakeTrigger.IsNull()) {
-    if (strcmp(gAlice->GetTriggerDescriptor(),""))
-      fMakeTrigger = gAlice->GetTriggerDescriptor();
-  }
-  else
-    gAlice->SetTriggerDescriptor(fMakeTrigger.Data());
   AliSysInfo::AddStamp("RunSimulation_TriggerDescriptor");
 
   // Set run number in CDBManager
@@ -1108,8 +1143,7 @@ Bool_t AliSimulation::RunSimulation(Int_t nEvents)
     loaderName += "Loader";
     AliLoader* loader = runLoader->GetLoader(loaderName);
     if (!loader) {
-      AliError(Form("RunSimulation", "no loader for %s found\n"
-                    "Number of events per file not set for %s %s", 
+      AliError(Form("RunSimulation no loader for %s found\n Number of events per file not set for %s %s", 
                     detName, typeName, detName));
       continue;
     }
@@ -1148,7 +1182,7 @@ Bool_t AliSimulation::RunSDigitization(const char* detectors)
 {
 // run the digitization and produce summable digits
   static Int_t eventNr=0;
-  AliCodeTimerAuto("") ;
+  AliCodeTimerAuto("",0) ;
 
   // initialize CDB storage, run number, set CDB lock
   InitCDB();
@@ -1190,7 +1224,7 @@ Bool_t AliSimulation::RunDigitization(const char* detectors,
 {
 // run the digitization and produce digits from sdigits
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
 
   // initialize CDB storage, run number, set CDB lock
   InitCDB();
@@ -1256,7 +1290,7 @@ Bool_t AliSimulation::RunHitsDigitization(const char* detectors)
 {
 // run the digitization and produce digits from hits
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
 
   // initialize CDB storage, run number, set CDB lock
   InitCDB();
@@ -1303,7 +1337,7 @@ Bool_t AliSimulation::WriteRawData(const char* detectors,
 // 'selrawdata' flag can be used to enable writing of detectors raw data
 // accoring to the trigger cluster.
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
   AliSysInfo::AddStamp("WriteRawData_Start");
   
   TString detStr = detectors;
@@ -1348,6 +1382,7 @@ Bool_t AliSimulation::WriteRawData(const char* detectors,
        sprintf(command, "rm -r raw%d", iEvent);
        gSystem->Exec(command);
       }
+      delete runLoader;
     }
 
     if (rootOutput) {
@@ -1379,7 +1414,7 @@ Bool_t AliSimulation::WriteRawFiles(const char* detectors)
 {
 // convert the digits to raw data DDL files
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
   
   AliRunLoader* runLoader = LoadRun("READ");
   if (!runLoader) return kFALSE;
@@ -1440,7 +1475,7 @@ Bool_t AliSimulation::ConvertRawFilesToDate(const char* dateFileName,
 // The second argument is not empty when the user decides to write
 // the detectors raw data according to the trigger cluster.
 
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
   
   char* path = gSystem->Which(gSystem->Getenv("PATH"), "dateStream");
   if (!path) {
@@ -1757,7 +1792,7 @@ Bool_t AliSimulation::IsSelected(TString detName, TString& detectors) const
 }
 
 //_____________________________________________________________________________
-Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* esdFileName) 
+Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* esdFileName, Int_t N
 {
 //
 // Steering routine  to convert raw data in directory rawDirectory/ to fake SDigits. 
@@ -1792,16 +1827,17 @@ Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* e
     
    AliPDG::AddParticlesToPdgDataBase();  
 
-   gAlice->GetMCApp()->Init();
-
    gMC->SetMagField(TGeoGlobalMagField::Instance()->GetField());
    
+   gAlice->GetMCApp()->Init();
+
    //Must be here because some MCs (G4) adds detectors here and not in Config.C
    gAlice->InitLoaders();
    AliRunLoader::Instance()->MakeTree("E");
    AliRunLoader::Instance()->LoadKinematics("RECREATE");
    AliRunLoader::Instance()->LoadTrackRefs("RECREATE");
    AliRunLoader::Instance()->LoadHits("all","RECREATE");
+
    //
    // Save stuff at the beginning of the file to avoid file corruption
    AliRunLoader::Instance()->CdGAFile();
@@ -1818,11 +1854,17 @@ Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* e
     AliRunLoader* runLoader = AliRunLoader::Instance();
     //
     // Open esd file if available
-    TFile* esdFile = TFile::Open(esdFileName);
+    TFile* esdFile = 0;
     TTree* treeESD = 0;
-    AliESDEvent* esd = new AliESDEvent();
-    esdFile->GetObject("esdTree", treeESD);
-    if (treeESD) esd->ReadFromTree(treeESD);
+    AliESDEvent* esd = 0;
+    if (esdFileName && (strlen(esdFileName)>0)) {
+      esdFile = TFile::Open(esdFileName);
+      if (esdFile) {
+        esd = new AliESDEvent();
+        esdFile->GetObject("esdTree", treeESD);
+        if (treeESD) esd->ReadFromTree(treeESD);
+      }
+    }
 
     //
     // Create the RawReader
@@ -1843,24 +1885,29 @@ Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* e
     //
     // Get Header
     AliHeader* header = runLoader->GetHeader();
-    //
-    TString detStr = fMakeSDigits;
     // Event loop
     Int_t nev = 0;
     while(kTRUE) {
        if (!(rawReader->NextEvent())) break;
+       runLoader->SetEventNumber(nev);
+        runLoader->GetHeader()->Reset(rawReader->GetRunNumber(), 
+                                      nev, nev);
+       runLoader->TreeE()->Fill();
+        runLoader->GetEvent(nev);
+        AliInfo(Form("We are at event %d",nev));
        //
        // Detector loop
+        TString detStr = fMakeSDigits;
        for (iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
            AliModule* det = (AliModule*) detArray->At(iDet);
            if (!det || !det->IsActive()) continue;
            if (IsSelected(det->GetName(), detStr)) {
-             AliInfo(Form("Calling Raw2SDigits for %s\n", det->GetName()));
+             AliInfo(Form("Calling Raw2SDigits for %s", det->GetName()));
              det->Raw2SDigits(rawReader);
              rawReader->Reset();
            }
        } // detectors
-
+        
 
        //
        //  If ESD information available obtain reconstructed vertex and store in header.
@@ -1876,13 +1923,14 @@ Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* e
            mcHeader->SetPrimaryVertex(mcV);
            header->Reset(0,nev);
            header->SetGenEventHeader(mcHeader);
-           printf("***** Saved vertex %f %f %f \n", position[0], position[1], position[2]);
+           AliInfo(Form("***** Saved vertex %f %f %f \n", position[0], position[1], position[2]));
        }
-       nev++;
 //
 //      Finish the event
-       runLoader->TreeE()->Fill();
-       runLoader->SetNextEvent();
+        AliInfo(Form("Finished event %d",nev));
+       nev++;
+        if (N>0&&nev>=N)
+          break;
     } // events
  
     delete rawReader;
@@ -2090,7 +2138,7 @@ Bool_t AliSimulation::RunQA()
 {
        // run the QA on summable hits, digits or digits
        
-  if(!gAlice) return kFALSE;
+    //if(!gAlice) return kFALSE;
        AliQAManager::QAManager()->SetRunLoader(AliRunLoader::Instance()) ;
 
        TString detectorsw("") ;  
@@ -2116,9 +2164,21 @@ Bool_t AliSimulation::SetRunQA(TString detAndAction)
        }
        Int_t colon = detAndAction.Index(":") ; 
        fQADetectors = detAndAction(0, colon) ; 
-       if (fQADetectors.Contains("ALL") )
-               fQADetectors = Form("%s %s", fMakeDigits.Data(), fMakeDigitsFromHits.Data()) ; 
-               fQATasks   = detAndAction(colon+1, detAndAction.Sizeof() ) ; 
+       if (fQADetectors.Contains("ALL") ){
+    TString tmp = Form("%s %s", fMakeDigits.Data(), fMakeDigitsFromHits.Data()) ; 
+    Int_t minus = fQADetectors.Last('-') ; 
+    TString toKeep = Form("%s %s", fMakeDigits.Data(), fMakeDigitsFromHits.Data()) ; 
+    TString toRemove("") ;
+    while (minus >= 0) {
+      toRemove = fQADetectors(minus+1, fQADetectors.Length()) ; 
+      toRemove = toRemove.Strip() ; 
+      toKeep.ReplaceAll(toRemove, "") ; 
+      fQADetectors.ReplaceAll(Form("-%s", toRemove.Data()), "") ; 
+      minus = fQADetectors.Last('-') ; 
+    }
+    fQADetectors = toKeep ; 
+  }
+  fQATasks   = detAndAction(colon+1, detAndAction.Sizeof() ) ; 
        if (fQATasks.Contains("ALL") ) {
                fQATasks = Form("%d %d %d", AliQAv1::kHITS, AliQAv1::kSDIGITS, AliQAv1::kDIGITS) ; 
        } else {
@@ -2204,13 +2264,17 @@ void AliSimulation::WriteGRPEntry()
 
   grpObj->SetRunType("PHYSICS");
   grpObj->SetTimeStart(0);
-  grpObj->SetTimeEnd(9999);
+  TDatime curtime;
+  grpObj->SetTimeStart(0);
+  grpObj->SetTimeEnd(curtime.Convert()); 
+  grpObj->SetBeamEnergyIsSqrtSHalfGeV(); // new format of GRP: store sqrt(s)/2 in GeV
 
   const AliGenerator *gen = gAlice->GetMCApp()->Generator();
+  Int_t a = 0;
+  Int_t z = 0;
+
   if (gen) {
-    grpObj->SetBeamEnergy(gen->GetEnergyCMS()/0.120);
     TString projectile;
-    Int_t a,z;
     gen->GetProjectile(projectile,a,z);
     TString target;
     gen->GetTarget(target,a,z);
@@ -2218,9 +2282,15 @@ void AliSimulation::WriteGRPEntry()
     beamType.ReplaceAll(" ","");
     if (!beamType.CompareTo("-")) {
       grpObj->SetBeamType("UNKNOWN");
+      grpObj->SetBeamEnergy(gen->GetEnergyCMS()/2);
     }
     else {
       grpObj->SetBeamType(beamType);
+      if (z != 0) {
+         grpObj->SetBeamEnergy(gen->GetEnergyCMS()/2 * a / z);
+      } else {
+         grpObj->SetBeamEnergy(gen->GetEnergyCMS()/2 );
+      }
       // Heavy ion run, the event specie is set to kHighMult
       fEventSpecie = AliRecoParam::kHighMult;
       if ((strcmp(beamType,"p-p") == 0) ||
@@ -2249,7 +2319,7 @@ void AliSimulation::WriteGRPEntry()
     }
   }
   // CTP
-  if (!fMakeTrigger.IsNull() || strcmp(gAlice->GetTriggerDescriptor(),""))
+  if (!fTriggerConfig.IsNull())
     detectorPattern |= (1 << AliDAQ::DetectorID("TRG"));
 
   // HLT
@@ -2260,8 +2330,6 @@ void AliSimulation::WriteGRPEntry()
   grpObj->SetDetectorMask((Int_t)detectorPattern);
   grpObj->SetLHCPeriod("LHC08c");
   grpObj->SetLHCState("STABLE_BEAMS");
-  grpObj->SetLHCLuminosity(0,(AliGRPObject::Stats)0);
-  grpObj->SetBeamIntensity(0,(AliGRPObject::Stats)0);
   //
   AliMagF *field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
   Float_t solenoidField = field ? TMath::Abs(field->SolenoidField()) : 0;
@@ -2286,6 +2354,11 @@ void AliSimulation::WriteGRPEntry()
 
   // Now store the entry in OCDB
   AliCDBManager* man = AliCDBManager::Instance();
+  
+  man->SetLock(0, fKey);
+  
+  AliCDBStorage* sto = man->GetStorage(fGRPWriteLocation.Data());
+  
 
   AliCDBId id("GRP/GRP/Data", man->GetRun(), man->GetRun(), 1, 1);
   AliCDBMetaData *metadata= new AliCDBMetaData();
@@ -2293,7 +2366,8 @@ void AliSimulation::WriteGRPEntry()
   metadata->SetResponsible("alice-off@cern.ch");
   metadata->SetComment("Automatically produced GRP entry for Monte Carlo");
  
-  man->Put(grpObj,id,metadata);
+  sto->Put(grpObj,id,metadata);
+  man->SetLock(1, fKey);
 }