[u/mrichter/AliRoot.git] / MUON / AliMUONReAlignTask.cxx

/**************************************************************************
* 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.                  *
**************************************************************************/

// $Id$

//-----------------------------------------------------------------------------
/// \class AliMUONReAlignTask
/// AliAnalysisTask to realign the MUON spectrometer.
/// The Task reads as input ESDs moves the clusters of a MUONTrack acoording
/// to the re aligned geometry taken from a misalignment file in the OCDB 
/// and refit the track. Then it produces a AliMUONClusterInfo object for each
/// cluster. The output is a TTree of AliMUONClusterInfo
///
/// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
//-----------------------------------------------------------------------------

#include <fstream>

#include <TString.h>
#include <TError.h>
#include <TTree.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TRandom.h>
#include <TGeoGlobalMagField.h>
#include <TGeoManager.h>
#include <Riostream.h>

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliESDInputHandler.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliMagF.h"
#include "AliCDBManager.h"
#include "AliGRPManager.h"
#include "AliGeomManager.h"

#include "AliMpConstants.h"
#include "AliMpCDB.h"
#include "AliMpSegmentation.h"
#include "AliMpVSegmentation.h"
#include "AliMpPad.h"
#include "AliMUONCalibrationData.h"
#include "AliMUONVCalibParam.h"
#include "AliMUONPadInfo.h"
#include "AliMUONClusterInfo.h"
#include "AliMUONRecoParam.h"
#include "AliMUONESDInterface.h"
#include "AliMUONRefitter.h"
#include "AliMUONRecoParam.h"
#include "AliMUONVDigit.h"
#include "AliMUONVCluster.h"
#include "AliMUONTrack.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONVDigitStore.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTrackExtrap.h"
#include "AliMUONGeometryTransformer.h"

#include "AliMUONReAlignTask.h"

///\cond CLASSIMP   
ClassImp(AliMUONReAlignTask)
///\endcond

//________________________________________________________________________
AliMUONReAlignTask::AliMUONReAlignTask(const char *name, const char *geofilename, const char *defaultocdb, const char *misalignocdb) 
  : AliAnalysisTask(name, ""),
    fESD(0x0),
    fClusterInfoTree(0x0),
    fClusterInfo(0x0),
    fESDInterface(0x0),
    fRefitter(0x0),
    fRecoParam(0x0),
    fGeoFilename(geofilename),
    fMisAlignOCDB(misalignocdb),
    fDefaultOCDB(defaultocdb),
    fGeoTransformer(0x0),
    fNewGeoTransformer(0x0),
    fGainStore(0x0),
    fPedStore(0x0),
    fPrintLevel(0),
    fLastRun(-1)
{
  /// Default Constructor
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes a TTree
  DefineOutput(0, TTree::Class());

  fClusterInfo = new AliMUONClusterInfo();
  fESDInterface = new AliMUONESDInterface();
  fGeoTransformer = new AliMUONGeometryTransformer();
  fNewGeoTransformer = new AliMUONGeometryTransformer();
}

//________________________________________________________________________
AliMUONReAlignTask::AliMUONReAlignTask(const AliMUONReAlignTask& obj)
  : AliAnalysisTask(obj),
    fESD(0x0),
    fClusterInfoTree(0x0),
    fClusterInfo(0x0),
    fESDInterface(0x0),
    fRefitter(0x0),
    fRecoParam(0x0),
    fGeoFilename(""),
    fMisAlignOCDB(""),
    fDefaultOCDB(""),
    fGeoTransformer(0x0),
    fNewGeoTransformer(0x0),
    fGainStore(0x0),
    fPedStore(0x0),
    fPrintLevel(0),
    fLastRun(-1)
{
  /// Copy constructor
  fESD = obj.fESD;
  fClusterInfoTree = obj.fClusterInfoTree;
  fClusterInfo = obj.fClusterInfo;
  fESDInterface = obj.fESDInterface;
  fRefitter = obj.fRefitter;
  fRecoParam = obj.fRecoParam;
  fGeoFilename = obj.fGeoFilename;
  fMisAlignOCDB = obj.fMisAlignOCDB;
  fDefaultOCDB = obj.fDefaultOCDB;
  fGeoTransformer = obj.fGeoTransformer;
  fNewGeoTransformer = obj.fNewGeoTransformer;
  fGainStore = obj.fGainStore;
  fPedStore = obj.fPedStore;
  fPrintLevel = obj.fPrintLevel;
  fLastRun = obj.fLastRun;
}

//________________________________________________________________________________
AliMUONReAlignTask& AliMUONReAlignTask::operator=(const AliMUONReAlignTask& other)
{
  /// Assignment
  if(&other == this) return *this;
  AliAnalysisTask::operator=(other);
  fESD = other.fESD;
  fClusterInfoTree = other.fClusterInfoTree;
  fClusterInfo = other.fClusterInfo;
  fESDInterface = other.fESDInterface;
  fRefitter = other.fRefitter;
  fRecoParam = other.fRecoParam;
  fGeoFilename = other.fGeoFilename;
  fMisAlignOCDB = other.fMisAlignOCDB;
  fDefaultOCDB = other.fDefaultOCDB;
  fGeoTransformer = other.fGeoTransformer;
  fNewGeoTransformer = other.fNewGeoTransformer;
  fGainStore = other.fGainStore;
  fPedStore = other.fPedStore;
  fPrintLevel = other.fPrintLevel;
  fLastRun = other.fLastRun;

  return *this;
} 


//________________________________________________________________________
AliMUONReAlignTask::~AliMUONReAlignTask() 
{ 
  /// Destructor
  if (fESDInterface) delete fESDInterface;
  if (fGeoTransformer) delete fGeoTransformer;
  if (fNewGeoTransformer) delete fNewGeoTransformer;
}

//________________________________________________________________________
void AliMUONReAlignTask::LocalInit() 
{
  /// Local initialization, called once per task on the client machine 
  /// where the analysis train is assembled
  AliMpCDB::LoadMpSegmentation();

  // prepare the refitting
  gRandom->SetSeed(0);
  Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());

  fRefitter = new AliMUONRefitter(fRecoParam);
  fRefitter->Connect(fESDInterface);
  
  // Original geotransformer
  fGeoTransformer->LoadGeometryData();   
  // Apply mis alignments
  AliGeomManager::ApplyAlignObjsFromCDB("MUON");    
  // New geotransformer
  fNewGeoTransformer->LoadGeometryData();   
  
}

//________________________________________________________________________
void AliMUONReAlignTask::ConnectInputData(Option_t *) 
{
  /// Connect ESD here. Called on each input data change.
  // Connect ESD here
  TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
  if (!esdTree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } 
  else {
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());   
    if (!esdH) {
      Printf("ERROR: Could not get ESDInputHandler");
    } 
    else {
      fESD = esdH->GetEvent();
    }
  }
}

//________________________________________________________________________
void AliMUONReAlignTask::CreateOutputObjects()
{
  /// Executed once on each worker (machine actually running the analysis code)
  //
  // This method has to be called INSIDE the user redefined CreateOutputObjects
  // method, before creating each object corresponding to the output containers
  // that are to be written to a file. This need to be done in general for the big output
  // objects that may not fit memory during processing. 
  OpenFile(0); 
  
  fClusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
  fClusterInfoTree->Branch("clusterInfo", &fClusterInfo, 32000, 99);

}

//________________________________________________________________________
void AliMUONReAlignTask::Exec(Option_t *) 
{
  /// Main loop, called for each event
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }

  // Prepare Gain and Pedestal store
  if (!fGainStore || fLastRun!=fESD->GetRunNumber()){
    fGainStore = AliMUONCalibrationData::CreateGains(fESD->GetRunNumber());
    fLastRun = fESD->GetRunNumber();
  }
  if (!fPedStore || fLastRun!=fESD->GetRunNumber()){
    fPedStore = AliMUONCalibrationData::CreatePedestals(fESD->GetRunNumber());
    fLastRun = fESD->GetRunNumber();
  }

  AliMUONPadInfo padInfo;

  Int_t nTracks = (Int_t)fESD->GetNumberOfMuonTracks();
  if (nTracks < 1) return;
  
  // load the current event
  fESDInterface->LoadEvent(*fESD);
  AliMUONVDigitStore* digitStore = fESDInterface->GetDigits();

  Double_t lX = 0.;
  Double_t lY = 0.;
  Double_t lZ = 0.;
  Double_t gX = 0.;
  Double_t gY = 0.;
  Double_t gZ = 0.;
  // loop over cluster to modify their position
  AliMUONVCluster *cluster;
  TIter next(fESDInterface->CreateClusterIterator());
  while ((cluster = static_cast<AliMUONVCluster*>(next()))) {
    //       cout << "Original cluster" << endl;
    //       cluster->Print();
    gX = cluster->GetX();
    gY = cluster->GetY();
    gZ = cluster->GetZ();
    fGeoTransformer->Global2Local(cluster->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
    fNewGeoTransformer->Local2Global(cluster->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
    cluster->SetXYZ(gX,gY,gZ);
    //       cout << "Aligned cluster" << endl;
    //       cluster->Print();
  }

  // refit the tracks from digits
  AliMUONVTrackStore* newTrackStore = fRefitter->ReconstructFromClusters();
    
  //----------------------------------------------//
  // ------ fill new ESD and print results ------ //
  //----------------------------------------------//
  // loop over the list of ESD tracks
  TClonesArray *esdTracks = (TClonesArray*) fESD->FindListObject("MuonTracks");
  for (Int_t iTrack = 0; iTrack <  nTracks; iTrack++) {      
    // get the ESD track
    AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
    
    // skip ghost tracks (leave them unchanged in the new ESD file)
    if (!esdTrack->ContainTrackerData()) continue;
      
    // get the corresponding MUON track
    AliMUONTrack* track = fESDInterface->FindTrack(esdTrack->GetUniqueID());
      
    // Find the corresponding re-fitted MUON track
    AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
    
//     // replace the content of the current ESD track or remove it if the refitting has failed
//     if (newTrack) {
//       Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
//       AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
//       } else {
//       esdTracks->Remove(esdTrack);
//       }
      
    // print initial and re-fitted track parameters at first cluster if any
    if (fPrintLevel>0) {
      cout<<"            ----------------track #"<<iTrack+1<<"----------------"<<endl;
      cout<<"before refit:"<<endl;
      AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
      param->Print("FULL");
      if (fPrintLevel>1) param->GetCovariances().Print();
      if (!newTrack) continue;
      cout<<"after refit:"<<endl;
      param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
      param->Print("FULL");
      if (fPrintLevel>1) param->GetCovariances().Print();
      cout<<"            ----------------------------------------"<<endl;
    }
    
    // Cluster Info part
    UInt_t muonClusterMap = BuildClusterMap(*newTrack);
	
    // loop over clusters
    AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->First());
    while (trackParam) {
      fClusterInfo->Clear("C");
      
      // fill cluster info
      cluster = trackParam->GetClusterPtr();
      fClusterInfo->SetRunId(fESD->GetRunNumber());
      fClusterInfo->SetEventId(fESD->GetEventNumberInFile());
      fClusterInfo->SetZ(cluster->GetZ());
      fClusterInfo->SetClusterId(cluster->GetUniqueID());
      fClusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
      fClusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
      fClusterInfo->SetClusterChi2(cluster->GetChi2());
      fClusterInfo->SetClusterCharge(cluster->GetCharge());
      
      // fill track info
      fClusterInfo->SetTrackId(newTrack->GetUniqueID());
      fClusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
      fClusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
      fClusterInfo->SetTrackP(trackParam->P());
      const TMatrixD paramCov = trackParam->GetCovariances();
      fClusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
      fClusterInfo->SetTrackChi2(newTrack->GetNormalizedChi2());
      fClusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
      fClusterInfo->SetTrackNHits(newTrack->GetNClusters());
      fClusterInfo->SetTrackChamberHitMap(muonClusterMap);
      
      // fill pad info if available	  
      for (Int_t i=0; i<cluster->GetNDigits(); i++) {
	AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
	if (!digit) continue;
	
	// pad location
	const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
	  GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
	AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
	
	// calibration parameters
	AliMUONVCalibParam* ped =  fPedStore ? static_cast<AliMUONVCalibParam*>(fPedStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
	AliMUONVCalibParam* gain = fGainStore ? static_cast<AliMUONVCalibParam*>(fGainStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
	Int_t manuChannel = digit->ManuChannel();
	Int_t planeType = 0;
	if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
	  planeType = 1;
	}
	
	// fill pad info
	padInfo.SetPadId(digit->GetUniqueID());
	padInfo.SetPadPlaneType(planeType);
	padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
	padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
	padInfo.SetPadCharge((Double_t)digit->Charge());
	padInfo.SetPadADC(digit->ADC());
	padInfo.SetSaturated(digit->IsSaturated());
	padInfo.SetCalibrated(digit->IsCalibrated());
	if (ped) {
	  padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
	} else {
	  padInfo.SetPedestal(-250.,-5.);
	}
	if (gain) {
	  padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
			  gain->ValueAsIntFast(manuChannel,2), gain->ValueAsIntFast(manuChannel,3));
	} else {
	  padInfo.SetGain(-1.,-0.1,-4095,-1);
	} 	
	fClusterInfo->AddPad(padInfo);
      }
      	  
      // fill cluster info tree
      fClusterInfoTree->Fill();
      trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
    }
  }
  // free memory
  delete newTrackStore;
    
  
  // Post final data. Write histo list to a file with option "RECREATE"
  PostData(0,fClusterInfoTree);
}      

//________________________________________________________________________
void AliMUONReAlignTask::Terminate(const Option_t*)
{
  /// Called once per task on the client machine at the end of the analysis.

}

//-----------------------------------------------------------------------
void AliMUONReAlignTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
{
  /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
  
  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
  if (!gGeoManager) {
    AliGeomManager::LoadGeometry(geoFilename);
    if (!gGeoManager) {
      Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
      return;
    }
  }
    
  // Load mapping
  AliCDBManager* man = AliCDBManager::Instance();
  man->SetDefaultStorage(defaultOCDB);
  man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
  man->Print();
  man->SetRun(0);
  if ( ! AliMpCDB::LoadDDLStore() ) {
    Error("MUONRefit","Could not access mapping from OCDB !");
    exit(-1);
  }

  // set mag field
  if (!TGeoGlobalMagField::Instance()->GetField()) {
    printf("Loading field map...\n");
    AliGRPManager *grpMan = new AliGRPManager();
    grpMan->ReadGRPEntry();
    grpMan->SetMagField();
    delete grpMan;
  }
  // set the magnetic field for track extrapolations
  AliMUONTrackExtrap::SetField();
  
  // Load initial reconstruction parameters from OCDB
  // reconstruction parameters
  fRecoParam = AliMUONRecoParam::GetCosmicParam();
  
  // digit selection
  fRecoParam->SetPadGoodnessMask(0x400BE80);
//   TString caliboption = caliboption1;
//   if ( calib == 2 ) caliboption = caliboption2;
//   fRecoParam->SetCalibrationMode("NOGAIN"caliboption.Data());
  fRecoParam->SetCalibrationMode("NOGAIN");
  
  // chamber resolution (incuding misalignment)
  for (Int_t iCh=0; iCh<10; iCh++) {
    fRecoParam->SetDefaultNonBendingReso(iCh,0.4);
    fRecoParam->SetDefaultBendingReso(iCh,0.4);
  }
  fRecoParam->SetMaxNonBendingDistanceToTrack(10.);
  fRecoParam->SetMaxBendingDistanceToTrack(10.);
  
  // cut on (non)bending slopes
  //fRecoParam->SetMaxNonBendingSlope(0.6);
  //fRecoParam->SetMaxBendingSlope(0.6);
  
  // tracking algorithm
  //  fRecoParam->MakeMoreTrackCandidates(kTRUE);
  fRecoParam->RequestStation(0, kFALSE);
  fRecoParam->RequestStation(2, kFALSE);
//   fRecoParam->RequestStation(3, kFALSE);
//   fRecoParam->RequestStation(4, kFALSE);
  fRecoParam->SetSigmaCutForTracking(7.);
  fRecoParam->ImproveTracks(kTRUE, 7.);
  Info("MUONRefit", "\n initial recontruction parameters:");
  fRecoParam->Print("FULL");
  
  AliMUONESDInterface::ResetTracker(fRecoParam);  
}

//-----------------------------------------------------------------------
UInt_t AliMUONReAlignTask::BuildClusterMap(AliMUONTrack &track)
{
  /// Build the map of clusters in tracking chambers
  
  UInt_t muonClusterMap = 0;
  
  AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
  while (trackParam) {
    
    muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
    
    trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
  }
  
  return muonClusterMap;
  
}
Commit	Line	Data
81f1d3ae	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// $Id$
	17
4d610fd5	18	//-----------------------------------------------------------------------------
	19	/// \class AliMUONReAlignTask
	20	/// AliAnalysisTask to realign the MUON spectrometer.
	21	/// The Task reads as input ESDs moves the clusters of a MUONTrack acoording
	22	/// to the re aligned geometry taken from a misalignment file in the OCDB
	23	/// and refit the track. Then it produces a AliMUONClusterInfo object for each
	24	/// cluster. The output is a TTree of AliMUONClusterInfo
	25	///
	26	/// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
	27	//-----------------------------------------------------------------------------
	28
	29	#include <fstream>
	30
	31	#include <TString.h>
	32	#include <TError.h>
	33	#include <TTree.h>
	34	#include <TChain.h>
	35	#include <TClonesArray.h>
	36	#include <TRandom.h>
	37	#include <TGeoGlobalMagField.h>
	38	#include <TGeoManager.h>
	39	#include <Riostream.h>
	40
	41	#include "AliAnalysisTask.h"
	42	#include "AliAnalysisManager.h"
	43	#include "AliESDInputHandler.h"
	44	#include "AliESDEvent.h"
	45	#include "AliESDMuonTrack.h"
	46	#include "AliMagF.h"
	47	#include "AliCDBManager.h"
cd8521dd	48	#include "AliGRPManager.h"
4d610fd5	49	#include "AliGeomManager.h"
	50
	51	#include "AliMpConstants.h"
	52	#include "AliMpCDB.h"
	53	#include "AliMpSegmentation.h"
	54	#include "AliMpVSegmentation.h"
	55	#include "AliMpPad.h"
	56	#include "AliMUONCalibrationData.h"
	57	#include "AliMUONVCalibParam.h"
	58	#include "AliMUONPadInfo.h"
	59	#include "AliMUONClusterInfo.h"
	60	#include "AliMUONRecoParam.h"
	61	#include "AliMUONESDInterface.h"
	62	#include "AliMUONRefitter.h"
	63	#include "AliMUONRecoParam.h"
	64	#include "AliMUONVDigit.h"
	65	#include "AliMUONVCluster.h"
	66	#include "AliMUONTrack.h"
	67	#include "AliMUONVTrackStore.h"
	68	#include "AliMUONVDigitStore.h"
	69	#include "AliMUONTrackParam.h"
	70	#include "AliMUONTrackExtrap.h"
	71	#include "AliMUONGeometryTransformer.h"
	72
	73	#include "AliMUONReAlignTask.h"
	74
	75	///\cond CLASSIMP
	76	ClassImp(AliMUONReAlignTask)
	77	///\endcond
	78
	79	//________________________________________________________________________
	80	AliMUONReAlignTask::AliMUONReAlignTask(const char name, const char geofilename, const char defaultocdb, const char misalignocdb)
	81	: AliAnalysisTask(name, ""),
	82	fESD(0x0),
	83	fClusterInfoTree(0x0),
	84	fClusterInfo(0x0),
	85	fESDInterface(0x0),
	86	fRefitter(0x0),
	87	fRecoParam(0x0),
	88	fGeoFilename(geofilename),
	89	fMisAlignOCDB(misalignocdb),
	90	fDefaultOCDB(defaultocdb),
	91	fGeoTransformer(0x0),
	92	fNewGeoTransformer(0x0),
	93	fGainStore(0x0),
	94	fPedStore(0x0),
	95	fPrintLevel(0),
	96	fLastRun(-1)
	97	{
	98	/// Default Constructor
	99	// Define input and output slots here
	100	// Input slot #0 works with a TChain
	101	DefineInput(0, TChain::Class());
	102	// Output slot #0 writes a TTree
	103	DefineOutput(0, TTree::Class());
	104
	105	fClusterInfo = new AliMUONClusterInfo();
	106	fESDInterface = new AliMUONESDInterface();
	107	fGeoTransformer = new AliMUONGeometryTransformer();
	108	fNewGeoTransformer = new AliMUONGeometryTransformer();
	109	}
	110
	111	//________________________________________________________________________
	112	AliMUONReAlignTask::AliMUONReAlignTask(const AliMUONReAlignTask& obj)
113	: AliAnalysisTask(obj),
114	fESD(0x0),
115	fClusterInfoTree(0x0),
116	fClusterInfo(0x0),
117	fESDInterface(0x0),
118	fRefitter(0x0),
119	fRecoParam(0x0),
120	fGeoFilename(""),
121	fMisAlignOCDB(""),
122	fDefaultOCDB(""),
123	fGeoTransformer(0x0),
124	fNewGeoTransformer(0x0),
125	fGainStore(0x0),
126	fPedStore(0x0),
127	fPrintLevel(0),
128	fLastRun(-1)
129	{
130	/// Copy constructor
131	fESD = obj.fESD;
132	fClusterInfoTree = obj.fClusterInfoTree;
133	fClusterInfo = obj.fClusterInfo;
134	fESDInterface = obj.fESDInterface;
135	fRefitter = obj.fRefitter;
136	fRecoParam = obj.fRecoParam;
137	fGeoFilename = obj.fGeoFilename;
138	fMisAlignOCDB = obj.fMisAlignOCDB;
139	fDefaultOCDB = obj.fDefaultOCDB;
140	fGeoTransformer = obj.fGeoTransformer;
141	fNewGeoTransformer = obj.fNewGeoTransformer;
142	fGainStore = obj.fGainStore;
143	fPedStore = obj.fPedStore;
144	fPrintLevel = obj.fPrintLevel;
145	fLastRun = obj.fLastRun;
146	}
147
148	//________________________________________________________________________________
149	AliMUONReAlignTask& AliMUONReAlignTask::operator=(const AliMUONReAlignTask& other)
150	{
151	/// Assignment
e99fb5c9	152	if(&other == this) return *this;
4d610fd5	153	AliAnalysisTask::operator=(other);
	154	fESD = other.fESD;
	155	fClusterInfoTree = other.fClusterInfoTree;
	156	fClusterInfo = other.fClusterInfo;
	157	fESDInterface = other.fESDInterface;
	158	fRefitter = other.fRefitter;
	159	fRecoParam = other.fRecoParam;
	160	fGeoFilename = other.fGeoFilename;
	161	fMisAlignOCDB = other.fMisAlignOCDB;
	162	fDefaultOCDB = other.fDefaultOCDB;
	163	fGeoTransformer = other.fGeoTransformer;
	164	fNewGeoTransformer = other.fNewGeoTransformer;
	165	fGainStore = other.fGainStore;
	166	fPedStore = other.fPedStore;
	167	fPrintLevel = other.fPrintLevel;
	168	fLastRun = other.fLastRun;
	169
	170	return *this;
	171	}
	172
	173
	174	//________________________________________________________________________
	175	AliMUONReAlignTask::~AliMUONReAlignTask()
	176	{
	177	/// Destructor
	178	if (fESDInterface) delete fESDInterface;
	179	if (fGeoTransformer) delete fGeoTransformer;
	180	if (fNewGeoTransformer) delete fNewGeoTransformer;
	181	}
	182
	183	//________________________________________________________________________
	184	void AliMUONReAlignTask::LocalInit()
	185	{
	186	/// Local initialization, called once per task on the client machine
	187	/// where the analysis train is assembled
	188	AliMpCDB::LoadMpSegmentation();
	189
	190	// prepare the refitting
	191	gRandom->SetSeed(0);
	192	Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
	193
	194	fRefitter = new AliMUONRefitter(fRecoParam);
	195	fRefitter->Connect(fESDInterface);
	196
	197	// Original geotransformer
	198	fGeoTransformer->LoadGeometryData();
	199	// Apply mis alignments
	200	AliGeomManager::ApplyAlignObjsFromCDB("MUON");
	201	// New geotransformer
	202	fNewGeoTransformer->LoadGeometryData();
	203
	204	}
	205
	206	//________________________________________________________________________
	207	void AliMUONReAlignTask::ConnectInputData(Option_t *)
	208	{
	209	/// Connect ESD here. Called on each input data change.
	210	// Connect ESD here
	211	TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
	212	if (!esdTree) {
	213	Printf("ERROR: Could not read chain from input slot 0");
	214	}
	215	else {
	216	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
217	if (!esdH) {
218	Printf("ERROR: Could not get ESDInputHandler");
219	}
220	else {
221	fESD = esdH->GetEvent();
222	}
223	}
224	}
225
226	//________________________________________________________________________
227	void AliMUONReAlignTask::CreateOutputObjects()
228	{
229	/// Executed once on each worker (machine actually running the analysis code)
230	//
231	// This method has to be called INSIDE the user redefined CreateOutputObjects
232	// method, before creating each object corresponding to the output containers
233	// that are to be written to a file. This need to be done in general for the big output
234	// objects that may not fit memory during processing.
235	OpenFile(0);
236
237	fClusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
238	fClusterInfoTree->Branch("clusterInfo", &fClusterInfo, 32000, 99);
239
240	}
241
242	//________________________________________________________________________
243	void AliMUONReAlignTask::Exec(Option_t *)
244	{
245	/// Main loop, called for each event
246	if (!fESD) {
247	Printf("ERROR: fESD not available");
248	return;
249	}
250
251	// Prepare Gain and Pedestal store
252	if (!fGainStore \|\| fLastRun!=fESD->GetRunNumber()){
253	fGainStore = AliMUONCalibrationData::CreateGains(fESD->GetRunNumber());
254	fLastRun = fESD->GetRunNumber();
255	}
256	if (!fPedStore \|\| fLastRun!=fESD->GetRunNumber()){
257	fPedStore = AliMUONCalibrationData::CreatePedestals(fESD->GetRunNumber());
258	fLastRun = fESD->GetRunNumber();
259	}
260
261	AliMUONPadInfo padInfo;
262
263	Int_t nTracks = (Int_t)fESD->GetNumberOfMuonTracks();
264	if (nTracks < 1) return;
265
266	// load the current event
267	fESDInterface->LoadEvent(*fESD);
268	AliMUONVDigitStore* digitStore = fESDInterface->GetDigits();
269
270	Double_t lX = 0.;
271	Double_t lY = 0.;
272	Double_t lZ = 0.;
273	Double_t gX = 0.;
274	Double_t gY = 0.;
275	Double_t gZ = 0.;
276	// loop over cluster to modify their position
277	AliMUONVCluster *cluster;
278	TIter next(fESDInterface->CreateClusterIterator());
279	while ((cluster = static_cast<AliMUONVCluster*>(next()))) {
280	// cout << "Original cluster" << endl;
281	// cluster->Print();
282	gX = cluster->GetX();
283	gY = cluster->GetY();
284	gZ = cluster->GetZ();
285	fGeoTransformer->Global2Local(cluster->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
286	fNewGeoTransformer->Local2Global(cluster->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
287	cluster->SetXYZ(gX,gY,gZ);
288	// cout << "Aligned cluster" << endl;
289	// cluster->Print();
290	}
291
292	// refit the tracks from digits
293	AliMUONVTrackStore* newTrackStore = fRefitter->ReconstructFromClusters();
294
295	//----------------------------------------------//
296	// ------ fill new ESD and print results ------ //
297	//----------------------------------------------//
298	// loop over the list of ESD tracks
299	TClonesArray esdTracks = (TClonesArray) fESD->FindListObject("MuonTracks");
300	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
301	// get the ESD track
302	AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
303
304	// skip ghost tracks (leave them unchanged in the new ESD file)
305	if (!esdTrack->ContainTrackerData()) continue;
306
307	// get the corresponding MUON track
308	AliMUONTrack* track = fESDInterface->FindTrack(esdTrack->GetUniqueID());
309
310	// Find the corresponding re-fitted MUON track
311	AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
312
313	// // replace the content of the current ESD track or remove it if the refitting has failed
314	// if (newTrack) {
315	// Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
316	// AliMUONESDInterface::MUONToESD(newTrack, esdTrack, vertex, esdInterface.GetDigits());
317	// } else {
318	// esdTracks->Remove(esdTrack);
319	// }
320
321	// print initial and re-fitted track parameters at first cluster if any
322	if (fPrintLevel>0) {
323	cout<<" ----------------track #"<<iTrack+1<<"----------------"<<endl;
324	cout<<"before refit:"<<endl;
325	AliMUONTrackParam param = (AliMUONTrackParam) track->GetTrackParamAtCluster()->First();
326	param->Print("FULL");
327	if (fPrintLevel>1) param->GetCovariances().Print();
328	if (!newTrack) continue;
329	cout<<"after refit:"<<endl;
330	param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
331	param->Print("FULL");
332	if (fPrintLevel>1) param->GetCovariances().Print();
333	cout<<" ----------------------------------------"<<endl;
334	}
335
336	// Cluster Info part
337	UInt_t muonClusterMap = BuildClusterMap(*newTrack);
338
339	// loop over clusters
340	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->First());
341	while (trackParam) {
342	fClusterInfo->Clear("C");
343
344	// fill cluster info
345	cluster = trackParam->GetClusterPtr();
346	fClusterInfo->SetRunId(fESD->GetRunNumber());
347	fClusterInfo->SetEventId(fESD->GetEventNumberInFile());
348	fClusterInfo->SetZ(cluster->GetZ());
349	fClusterInfo->SetClusterId(cluster->GetUniqueID());
350	fClusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
351	fClusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
352	fClusterInfo->SetClusterChi2(cluster->GetChi2());
353	fClusterInfo->SetClusterCharge(cluster->GetCharge());
354
355	// fill track info
356	fClusterInfo->SetTrackId(newTrack->GetUniqueID());
357	fClusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
358	fClusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
359	fClusterInfo->SetTrackP(trackParam->P());
360	const TMatrixD paramCov = trackParam->GetCovariances();
361	fClusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
362	fClusterInfo->SetTrackChi2(newTrack->GetNormalizedChi2());
363	fClusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
364	fClusterInfo->SetTrackNHits(newTrack->GetNClusters());
365	fClusterInfo->SetTrackChamberHitMap(muonClusterMap);
366
367	// fill pad info if available
368	for (Int_t i=0; i<cluster->GetNDigits(); i++) {
369	AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
370	if (!digit) continue;
371
372	// pad location
373	const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
374	GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
375	AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
376
377	// calibration parameters
378	AliMUONVCalibParam* ped = fPedStore ? static_cast<AliMUONVCalibParam*>(fPedStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
379	AliMUONVCalibParam* gain = fGainStore ? static_cast<AliMUONVCalibParam*>(fGainStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
380	Int_t manuChannel = digit->ManuChannel();
381	Int_t planeType = 0;
382	if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
383	planeType = 1;
384	}
385
386	// fill pad info
387	padInfo.SetPadId(digit->GetUniqueID());
388	padInfo.SetPadPlaneType(planeType);
389	padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
390	padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
391	padInfo.SetPadCharge((Double_t)digit->Charge());
392	padInfo.SetPadADC(digit->ADC());
393	padInfo.SetSaturated(digit->IsSaturated());
394	padInfo.SetCalibrated(digit->IsCalibrated());
f4c2989f	395	if (ped) {
	396	padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
	397	} else {
	398	padInfo.SetPedestal(-250.,-5.);
	399	}
	400	if (gain) {
	401	padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
	402	gain->ValueAsIntFast(manuChannel,2), gain->ValueAsIntFast(manuChannel,3));
	403	} else {
	404	padInfo.SetGain(-1.,-0.1,-4095,-1);
	405	}
4d610fd5	406	fClusterInfo->AddPad(padInfo);
	407	}
	408
	409	// fill cluster info tree
	410	fClusterInfoTree->Fill();
	411	trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
	412	}
	413	}
	414	// free memory
	415	delete newTrackStore;
	416
	417
	418	// Post final data. Write histo list to a file with option "RECREATE"
	419	PostData(0,fClusterInfoTree);
	420	}
	421
	422	//________________________________________________________________________
	423	void AliMUONReAlignTask::Terminate(const Option_t*)
	424	{
	425	/// Called once per task on the client machine at the end of the analysis.
	426
	427	}
	428
	429	//-----------------------------------------------------------------------
	430	void AliMUONReAlignTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
	431	{
	432	/// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
	433
	434	// Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
	435	if (!gGeoManager) {
	436	AliGeomManager::LoadGeometry(geoFilename);
	437	if (!gGeoManager) {
	438	Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
	439	return;
	440	}
	441	}
cd8521dd	442
4d610fd5	443	// Load mapping
	444	AliCDBManager* man = AliCDBManager::Instance();
	445	man->SetDefaultStorage(defaultOCDB);
	446	man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
	447	man->Print();
	448	man->SetRun(0);
	449	if ( ! AliMpCDB::LoadDDLStore() ) {
	450	Error("MUONRefit","Could not access mapping from OCDB !");
	451	exit(-1);
	452	}
cd8521dd	453
	454	// set mag field
	455	if (!TGeoGlobalMagField::Instance()->GetField()) {
	456	printf("Loading field map...\n");
	457	AliGRPManager *grpMan = new AliGRPManager();
	458	grpMan->ReadGRPEntry();
	459	grpMan->SetMagField();
	460	delete grpMan;
	461	}
	462	// set the magnetic field for track extrapolations
	463	AliMUONTrackExtrap::SetField();
4d610fd5	464
	465	// Load initial reconstruction parameters from OCDB
	466	// reconstruction parameters
	467	fRecoParam = AliMUONRecoParam::GetCosmicParam();
	468
	469	// digit selection
	470	fRecoParam->SetPadGoodnessMask(0x400BE80);
	471	// TString caliboption = caliboption1;
	472	// if ( calib == 2 ) caliboption = caliboption2;
	473	// fRecoParam->SetCalibrationMode("NOGAIN"caliboption.Data());
	474	fRecoParam->SetCalibrationMode("NOGAIN");
	475
	476	// chamber resolution (incuding misalignment)
	477	for (Int_t iCh=0; iCh<10; iCh++) {
	478	fRecoParam->SetDefaultNonBendingReso(iCh,0.4);
	479	fRecoParam->SetDefaultBendingReso(iCh,0.4);
	480	}
	481	fRecoParam->SetMaxNonBendingDistanceToTrack(10.);
	482	fRecoParam->SetMaxBendingDistanceToTrack(10.);
	483
	484	// cut on (non)bending slopes
	485	//fRecoParam->SetMaxNonBendingSlope(0.6);
	486	//fRecoParam->SetMaxBendingSlope(0.6);
	487
	488	// tracking algorithm
	489	// fRecoParam->MakeMoreTrackCandidates(kTRUE);
	490	fRecoParam->RequestStation(0, kFALSE);
	491	fRecoParam->RequestStation(2, kFALSE);
	492	// fRecoParam->RequestStation(3, kFALSE);
	493	// fRecoParam->RequestStation(4, kFALSE);
	494	fRecoParam->SetSigmaCutForTracking(7.);
	495	fRecoParam->ImproveTracks(kTRUE, 7.);
	496	Info("MUONRefit", "\n initial recontruction parameters:");
	497	fRecoParam->Print("FULL");
	498
	499	AliMUONESDInterface::ResetTracker(fRecoParam);
	500	}
	501
	502	//-----------------------------------------------------------------------
	503	UInt_t AliMUONReAlignTask::BuildClusterMap(AliMUONTrack &track)
	504	{
	505	/// Build the map of clusters in tracking chambers
	506
	507	UInt_t muonClusterMap = 0;
	508
	509	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
	510	while (trackParam) {
	511
	512	muonClusterMap \|= BIT(trackParam->GetClusterPtr()->GetChamberId());
	513
	514	trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
	515	}
	516
	517	return muonClusterMap;
	518
	519	}