]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/MUONAlignment.C
Histogram ranges changed to cut off saturation peak and noise
[u/mrichter/AliRoot.git] / MUON / MUONAlignment.C
index 26a18c9a978c037da62d9f6ac54c2f4cf8254619..9d7d636d7b75a5ee66002f2cb21f0db0b670d421 100644 (file)
 
 /* $Id$ */
 
-// ---
-// Macro for MUON alignment using physics tracks. The macro uses AliMUONAlignment
-// class to calculate the alignment parameters.
-// An array for the alignment parameters is created and can be filled with
-// initial values that will be used as starting values by the alignment
-// algorithm.
-// Ny default the macro run over galice.root in the working directory. If a file list
-// of galice.root is provided as third argument the macro will run over all files.
-// The macro loop over the files, events and tracks. For each track
-// AliMUONAlignment::ProcessTrack(AliMUONTrack * track) and then
-// AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t
-// lSingleFit) are called. After all tracks have been procesed and fitted
-// AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls)
-// is called. The array parameters contains the obatained misalignement parameters.
-// A realigned geometry is generated in a local CDB.
-//
-// Authors: B. Becker and J. Castillo
-// ---
-
-void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString fileName = "galice.root", TString fileList = "")
+/// \ingroup macros
+/// \file MUONAlignment.C
+/// \brief Macro for MUON alignment using physics tracks. 
+///
+/// The macro uses the AliMUONAlignment class to calculate the alignment parameters.
+/// An array for the alignment parameters is created and can be filled with
+/// initial values that will be used as starting values by the alignment
+/// algorithm.
+///
+/// By default the macro run over galice.root in the working directory. If a file list
+/// of galice.root is provided as third argument the macro will run over all files.
+/// The macro loop over the files, events and tracks. For each track
+/// AliMUONAlignment::ProcessTrack(AliMUONTrack * track) and then
+/// AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t
+/// lSingleFit) are called. After all tracks have been procesed and fitted
+/// AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls)
+/// is called. The array parameters contains the obatained misalignement parameters.
+/// A realigned geometry is generated in a local CDB.
+///
+/// \author: B. Becker and J. Castillo
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "AliMUONAlignment.h"
+#include "AliMUONTrack.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONESDInterface.h"
+
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliMagFMaps.h"
+#include "AliTracker.h"
+#include "AliCDBManager.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+#include "AliGeomManager.h"
+
+#include <TString.h>
+#include <TError.h>
+#include <TH1.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TClonesArray.h>
+#include <Riostream.h>
+
+#include <fstream>
+
+#endif
+
+void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString esdFileName = "AliESDs.root", TString fileList = "")
 {
  
   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
-  if (!gGeoManager) {
-    TGeoManager::Import(geoFilename);
-    if (!gGeoManager) {
-      Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
+  if ( ! AliGeomManager::GetGeometry() ) {
+    AliGeomManager::LoadGeometry(geoFilename);
+    if (! AliGeomManager::GetGeometry() ) {
+      Error("MUONAlignment", "getting geometry from file %s failed", geoFilename);
       return;
     }
   }
@@ -65,10 +98,10 @@ void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root",
 
   Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
 
-  Int_t iPar = 0;
   // Set initial values here, good guess may help convergence
   // St 1 
-  //   parameters[iPar++] =  0.010300 ;  parameters[iPar++] =  0.010600 ;  parameters[iPar++] =  0.000396 ;  
+  //  Int_t iPar = 0;
+  //  parameters[iPar++] =  0.010300 ;  parameters[iPar++] =  0.010600 ;  parameters[iPar++] =  0.000396 ;  
 
   bool bLoop = kFALSE;
   ifstream sFileList;
@@ -88,16 +121,40 @@ void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root",
   AliMUONAlignment* alig = new AliMUONAlignment();
   alig->InitGlobalParameters(parameters);
 
-  AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer(true);
-  transform->ReadGeometryData("volpath.dat", gGeoManager);
+  AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
+  transform->LoadGeometryData();
   alig->SetGeometryTransformer(transform);
 
-  char cFileName[100];  
-  AliMUONDataInterface amdi;
+  // Do alignment with magnetic field off
+  alig->SetBFieldOn(kFALSE);
+  
+  // Set tracking station to use
+  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
+
+  // Fix parameters or add constraints here
+  for (Int_t iSt=0; iSt<5; iSt++)
+    if (!bStOnOff[iSt]) alig->FixStation(iSt+1);
+
+  // Left and right sides of the detector are independent, one can choose to align 
+  // only one side
+  Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
+  alig->FixHalfSpectrometer(bStOnOff,bSpecLROnOff);
+
+  // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
+  Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
+  Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
+  alig->AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
 
+
+  char cFileName[100];  
+    
   Int_t lMaxFile = 1000;
   Int_t iFile = 0;
+  Int_t iEvent = 0;
   bool bKeepLoop = kTRUE;
+  Int_t iTrackTot=0;
+  Int_t iTrackOk=0;
+
   while(bKeepLoop && iFile<lMaxFile){
     iFile++;
     if (bLoop) {
@@ -105,36 +162,61 @@ void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root",
       if (sFileList.eof()) bKeepLoop = kFALSE;
     }
     else {
-      sprintf(cFileName,fileName.Data());
+      sprintf(cFileName,esdFileName.Data());
       bKeepLoop = kFALSE;
     }
-    if (!strstr(cFileName,"galice.root")) continue;      
+    if (!strstr(cFileName,"AliESDs.root")) continue;      
     cout << "Using file: " << cFileName << endl;
-    amdi.SetFile(cFileName);
-    Int_t nevents = amdi.NumberOfEvents();
+    
+    // load ESD event
+    TFile* esdFile = TFile::Open(cFileName); // open the file
+    if (!esdFile || !esdFile->IsOpen()) {
+      cout << "opening ESD file " << cFileName << "failed" << endl;
+      continue;
+    }
+    TTree* esdTree = (TTree*) esdFile->Get("esdTree"); // get the tree
+    if (!esdTree) {
+      cout << "no ESD tree found" << endl;
+      esdFile->Close();
+      continue;
+    }
+    AliESDEvent* esdEvent = new AliESDEvent(); // link ESD event to the tree
+    esdEvent->ReadFromTree(esdTree);
+
+    Int_t nevents = esdTree->GetEntries();
+    cout << "... with " << nevents << endl;
     for(Int_t event = 0; event < nevents; event++) {
-      amdi.GetEvent(event);
-      Int_t ntracks = amdi.NumberOfRecTracks();
-      cout << " there are " << ntracks << " tracks in event " << event << endl;
-      Int_t iTrack=0;
-      Int_t iTrackOk=0;
-      AliMUONTrack* track = amdi.RecTrack(iTrack);
-      while(track) {   
-       AliMUONTrackParam trackParam(*((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First())));
-       AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.);
-       Double_t invBenMom = trackParam.GetInverseBendingMomentum();
+      if (iEvent >= nEvents){
+       bKeepLoop = kFALSE;
+       break;
+      }
+      iEvent++;
+
+      if (esdTree->GetEvent(event) <= 0) {
+       cout << "fails to read ESD object for event " << event << endl;
+       continue;
+      }
+
+      Int_t nTracks = Int_t(esdEvent->GetNumberOfMuonTracks());
+      if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
+      for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+       AliESDMuonTrack* esdTrack = esdEvent->GetMuonTrack(iTrack);
+       if (!esdTrack->ClustersStored()) continue;
+       Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
        fInvBenMom->Fill(invBenMom);
        fBenMom->Fill(1./invBenMom);
        if (TMath::Abs(invBenMom)<=1.04) {
-         cout << "Track " << iTrack << endl;
-         alig->ProcessTrack(track);
-         cout << "Calling LocalFit" << endl;
+         AliMUONTrack track;
+         AliMUONESDInterface::ESDToMUON(*esdTrack, track);
+         alig->ProcessTrack(&track);
          alig->LocalFit(iTrackOk++,trackParams,0);
        }
-       iTrack++;
-       track = amdi.RecTrack(iTrack);
+       iTrackTot++;
       }
     }
+    delete esdEvent;
+    esdFile->Close();
+    cout << "Processed " << iTrackTot << " Tracks so far." << endl;
   }
   alig->GlobalFit(parameters,errors,pulls);
 
@@ -153,7 +235,7 @@ void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root",
   Double_t MSDEytErr[156] = {0};
   Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
   Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
-  Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
+  // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
   Int_t idOffset = 0; // 400
   Int_t lSDetElemCh = 0;
   for(Int_t iDE=0; iDE<lNDetElem; iDE++){
@@ -205,7 +287,7 @@ void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root",
   newTransform->WriteTransformations("transform2ReAlign.dat");
   
   // Generate realigned data in local cdb
-  TClonesArray* array = newTransform->GetMisAlignmentData();
+  const TClonesArray* array = newTransform->GetMisAlignmentData();
    
   // CDB manager
   AliCDBManager* cdbManager = AliCDBManager::Instance();
@@ -215,6 +297,7 @@ void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root",
   cdbData->SetResponsible("Dimuon Offline project");
   cdbData->SetComment("MUON alignment objects with residual misalignment");
   AliCDBId id("MUON/Align/Data", 0, 0); 
-  cdbManager->Put(array, id, cdbData);
+  cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
+
+}
 
-}