[u/mrichter/AliRoot.git] / MUON / MUONClusterInfo.C

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

/// \ingroup macros
/// \file MUONClusterInfo.C
/// \brief Macro to fill AliMUONClusterInfo objects
///
/// \author Philippe Pillot, SUBATECH

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TStopwatch.h>
#include <TFile.h>
#include <TClonesArray.h>
#include <TTree.h>
#include <TString.h>
#include <Riostream.h>
#include <TGeoManager.h>
#include <TRandom.h>
#include <TROOT.h>

// STEER includes
#include "AliMagFMaps.h"
#include "AliTracker.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliRecoParam.h"
#include "AliCDBManager.h"
#include "AliGeomManager.h"

// MUON includes
#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 "AliMUONVDigit.h"
#include "AliMUONVDigitStore.h"
#include "AliMUONVCluster.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackParam.h"
#endif

const Int_t printLevel = 1;

void Prepare();
TTree* GetESDTree(TFile *esdFile);

//-----------------------------------------------------------------------
void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* outFileName = "clusterInfo.root")
{
  /// load ESD event in the ESDInterface to recover MUON objects;
  /// fill AliMUONESDClusterInfo object with ESD cluster + corresponding track parameters;
  /// write results in a new root file.
  
  AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
  AliMUONPadInfo padInfo;
  AliMUONCalibrationData* calibData = 0x0;
  AliMUONESDInterface esdInterface;
  
  // prepare the refitting during ESD->MUON conversion
  Prepare();
  
  // open the ESD file and tree and connect the ESD event
  TFile* esdFile = TFile::Open(esdFileName);
  TTree* esdTree = GetESDTree(esdFile);
  AliESDEvent* esd = new AliESDEvent();
  esd->ReadFromTree(esdTree);
  
  // prepare the output tree
  gROOT->cd();
  TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
  clusterInfoFile->SetCompressionLevel(1);
  
  TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
  clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
  
  // timer start...
  TStopwatch timer;
  
  // Loop over ESD events
  if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
  else nevents = (Int_t)esdTree->GetEntries();
  for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
    
    //----------------------------------------------//
    // -------------- process event --------------- //
    //----------------------------------------------//
    // get the ESD of current event
    esdTree->GetEvent(iEvent);
    if (!esd) {
      Error("CheckESD", "no ESD object found for event %d", iEvent);
      return;
    }
    if (esd->GetNumberOfMuonTracks() < 1) continue;
    
    // prepare access to calibration data
    if (!calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
    
    // load the current event
    esdInterface.LoadEvent(*esd);
    
    // get digit store
    AliMUONVDigitStore* digitStore = esdInterface.GetDigits();
    
    //----------------------------------------------//
    // ------------- fill cluster info ------------ //
    //----------------------------------------------//
    // loop over the refitted tracks
    TIter nextTrack(esdInterface.CreateTrackIterator());
    AliMUONTrack* track;
    while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
      
      // loop over clusters
      AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
      while (trackParam) {
	clusterInfo->Clear("C");
	
	// fill cluster info
	AliMUONVCluster* cluster = trackParam->GetClusterPtr();
	clusterInfo->SetEventId(iEvent);
	clusterInfo->SetZ(cluster->GetZ());
	clusterInfo->SetClusterId(cluster->GetUniqueID());
	clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
	clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
	clusterInfo->SetClusterChi2(cluster->GetChi2());
	clusterInfo->SetClusterCharge(cluster->GetCharge());
	
	// fill track info
	clusterInfo->SetTrackId(track->GetUniqueID());
	clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
	clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetBendingSlope()), TMath::ATan(trackParam->GetNonBendingSlope()));
	clusterInfo->SetTrackP(trackParam->P());
	const TMatrixD paramCov = trackParam->GetCovariances();
	clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
	clusterInfo->SetTrackChi2(track->GetNormalizedChi2());
	clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());

	// 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(AliMpIntPair(digit->PadX(), digit->PadY()));
	  
	  // calibration parameters
	  AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
	  AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
	  Int_t manuChannel = digit->ManuChannel();
	  
	  // fill pad info
	  padInfo.SetPadId(digit->GetUniqueID());
	  padInfo.SetPadXY(pad.Position().X(), pad.Position().Y());
	  padInfo.SetPadDimXY(pad.Dimensions().X(), pad.Dimensions().Y());
	  padInfo.SetPadCharge((Double_t)digit->Charge());
	  padInfo.SetPadADC(digit->ADC());
	  padInfo.SetSaturated(digit->IsSaturated());
	  padInfo.SetCalibrated(digit->IsCalibrated());
	  padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
	  padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
			  gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
	  
	  clusterInfo->AddPad(padInfo);
	}
	
	// fill cluster info tree
	clusterInfoTree->Fill();
	
	trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
      }
      
    }
    
  }
  
  // ...timer stop
  timer.Stop();
  printf("Writing Tree\n");
  // write output tree
  clusterInfoFile->cd();
  clusterInfoTree->Write();
  printf("Deleting Tree\n");
  delete clusterInfoTree;
  printf("Closing File\n");
  clusterInfoFile->Close();
  
  // free memory
  printf("Deleting calibData\n");
  delete calibData;
  printf("Deleting clusterInfo\n");
  delete clusterInfo;
  printf("Closing esdFile\n");
  esdFile->Close();
  printf("Deleting esd\n");
  delete esd;
  //  delete padInfo;
  cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
}

//-----------------------------------------------------------------------
void Prepare()
{
  /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
  
  gRandom->SetSeed(0);
  
  // Import TGeo geometry (needed for track extrapolation)
  if (!gGeoManager) {
    AliGeomManager::LoadGeometry("geometry.root");
    if (!gGeoManager) {
      Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
      return;
    }
  }
  
  // set mag field
  printf("Loading field map...\n");
  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
  AliTracker::SetFieldMap(field, kFALSE);
  
  // Load mapping
  AliCDBManager* man = AliCDBManager::Instance();
  man->SetDefaultStorage("local://$ALICE_ROOT");
  man->SetRun(0);
  if ( ! AliMpCDB::LoadDDLStore() ) {
    Error("MUONRefit","Could not access mapping from OCDB !");
    exit(-1);
  }
  
  // eventually set reconstruction parameters for refit (otherwise read from OCDB)
/*  AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
  muonRecoParam->Print("FULL");
  AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
*/  
}

//-----------------------------------------------------------------------
TTree* GetESDTree(TFile *esdFile)
{
  /// Check that the file is properly open
  /// Return pointer to the ESD Tree
  
  if (!esdFile || !esdFile->IsOpen()) {
    Error("GetESDTree", "opening ESD file failed");
    exit(-1);
  }
  
  TTree* tree = (TTree*) esdFile->Get("esdTree");
  if (!tree) {
    Error("GetESDTree", "no ESD tree found");
    exit(-1);
  }
  
  return tree;
  
}
Commit	Line	Data
fcefbac4	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
	18	/// \ingroup macros
	19	/// \file MUONClusterInfo.C
	20	/// \brief Macro to fill AliMUONClusterInfo objects
	21	///
	22	/// \author Philippe Pillot, SUBATECH
	23
	24	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	25	#include <TStopwatch.h>
	26	#include <TFile.h>
	27	#include <TClonesArray.h>
	28	#include <TTree.h>
	29	#include <TString.h>
	30	#include <Riostream.h>
	31	#include <TGeoManager.h>
	32	#include <TRandom.h>
	33	#include <TROOT.h>
	34
	35	// STEER includes
	36	#include "AliMagFMaps.h"
	37	#include "AliTracker.h"
	38	#include "AliESDEvent.h"
	39	#include "AliESDMuonTrack.h"
	40	#include "AliRecoParam.h"
	41	#include "AliCDBManager.h"
	42	#include "AliGeomManager.h"
	43
	44	// MUON includes
	45	#include "AliMpCDB.h"
	46	#include "AliMpSegmentation.h"
	47	#include "AliMpVSegmentation.h"
	48	#include "AliMpPad.h"
	49	#include "AliMUONCalibrationData.h"
	50	#include "AliMUONVCalibParam.h"
	51	#include "AliMUONPadInfo.h"
	52	#include "AliMUONClusterInfo.h"
	53	#include "AliMUONRecoParam.h"
	54	#include "AliMUONESDInterface.h"
fcefbac4	55	#include "AliMUONVDigit.h"
	56	#include "AliMUONVDigitStore.h"
	57	#include "AliMUONVCluster.h"
fcefbac4	58	#include "AliMUONTrack.h"
fcefbac4	59	#include "AliMUONTrackParam.h"
fcefbac4	60	#endif
	61
	62	const Int_t printLevel = 1;
	63
	64	void Prepare();
	65	TTree* GetESDTree(TFile *esdFile);
	66
	67	//-----------------------------------------------------------------------
	68	void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* outFileName = "clusterInfo.root")
	69	{
8468d8c6	70	/// load ESD event in the ESDInterface to recover MUON objects;
fcefbac4	71	/// fill AliMUONESDClusterInfo object with ESD cluster + corresponding track parameters;
	72	/// write results in a new root file.
	73
	74	AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
	75	AliMUONPadInfo padInfo;
	76	AliMUONCalibrationData* calibData = 0x0;
8468d8c6	77	AliMUONESDInterface esdInterface;
fcefbac4	78
8468d8c6	79	// prepare the refitting during ESD->MUON conversion
fcefbac4	80	Prepare();
fcefbac4	81
	82	// open the ESD file and tree and connect the ESD event
	83	TFile* esdFile = TFile::Open(esdFileName);
	84	TTree* esdTree = GetESDTree(esdFile);
	85	AliESDEvent* esd = new AliESDEvent();
	86	esd->ReadFromTree(esdTree);
	87
	88	// prepare the output tree
	89	gROOT->cd();
	90	TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
	91	clusterInfoFile->SetCompressionLevel(1);
	92
	93	TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
	94	clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
	95
	96	// timer start...
	97	TStopwatch timer;
	98
	99	// Loop over ESD events
	100	if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
	101	else nevents = (Int_t)esdTree->GetEntries();
	102	for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
	103
	104	//----------------------------------------------//
	105	// -------------- process event --------------- //
	106	//----------------------------------------------//
	107	// get the ESD of current event
	108	esdTree->GetEvent(iEvent);
	109	if (!esd) {
	110	Error("CheckESD", "no ESD object found for event %d", iEvent);
	111	return;
	112	}
8468d8c6	113	if (esd->GetNumberOfMuonTracks() < 1) continue;
fcefbac4	114
	115	// prepare access to calibration data
	116	if (!calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
	117
	118	// load the current event
	119	esdInterface.LoadEvent(*esd);
	120
	121	// get digit store
	122	AliMUONVDigitStore* digitStore = esdInterface.GetDigits();
	123
fcefbac4	124	//----------------------------------------------//
	125	// ------------- fill cluster info ------------ //
	126	//----------------------------------------------//
	127	// loop over the refitted tracks
8468d8c6	128	TIter nextTrack(esdInterface.CreateTrackIterator());
	129	AliMUONTrack* track;
	130	while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
fcefbac4	131
fcefbac4	132	// loop over clusters
8468d8c6	133	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
fcefbac4	134	while (trackParam) {
	135	clusterInfo->Clear("C");
	136
	137	// fill cluster info
	138	AliMUONVCluster* cluster = trackParam->GetClusterPtr();
	139	clusterInfo->SetEventId(iEvent);
	140	clusterInfo->SetZ(cluster->GetZ());
	141	clusterInfo->SetClusterId(cluster->GetUniqueID());
	142	clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
	143	clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
	144	clusterInfo->SetClusterChi2(cluster->GetChi2());
	145	clusterInfo->SetClusterCharge(cluster->GetCharge());
	146
	147	// fill track info
8468d8c6	148	clusterInfo->SetTrackId(track->GetUniqueID());
fcefbac4	149	clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
	150	clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetBendingSlope()), TMath::ATan(trackParam->GetNonBendingSlope()));
	151	clusterInfo->SetTrackP(trackParam->P());
	152	const TMatrixD paramCov = trackParam->GetCovariances();
	153	clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
8468d8c6	154	clusterInfo->SetTrackChi2(track->GetNormalizedChi2());
fcefbac4	155	clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
	156
	157	// fill pad info if available
	158	for (Int_t i=0; i<cluster->GetNDigits(); i++) {
	159	AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
	160	if (!digit) continue;
	161
	162	// pad location
	163	const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
	164	GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
	165	AliMpPad pad = seg->PadByIndices(AliMpIntPair(digit->PadX(), digit->PadY()));
	166
	167	// calibration parameters
	168	AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
	169	AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
	170	Int_t manuChannel = digit->ManuChannel();
	171
	172	// fill pad info
	173	padInfo.SetPadId(digit->GetUniqueID());
	174	padInfo.SetPadXY(pad.Position().X(), pad.Position().Y());
	175	padInfo.SetPadDimXY(pad.Dimensions().X(), pad.Dimensions().Y());
	176	padInfo.SetPadCharge((Double_t)digit->Charge());
	177	padInfo.SetPadADC(digit->ADC());
	178	padInfo.SetSaturated(digit->IsSaturated());
	179	padInfo.SetCalibrated(digit->IsCalibrated());
	180	padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
	181	padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
	182	gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
	183
	184	clusterInfo->AddPad(padInfo);
	185	}
	186
	187	// fill cluster info tree
	188	clusterInfoTree->Fill();
	189
8468d8c6	190	trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
fcefbac4	191	}
	192
	193	}
	194
fcefbac4	195	}
	196
	197	// ...timer stop
	198	timer.Stop();
	199	printf("Writing Tree\n");
	200	// write output tree
	201	clusterInfoFile->cd();
	202	clusterInfoTree->Write();
	203	printf("Deleting Tree\n");
	204	delete clusterInfoTree;
	205	printf("Closing File\n");
	206	clusterInfoFile->Close();
	207
	208	// free memory
	209	printf("Deleting calibData\n");
	210	delete calibData;
	211	printf("Deleting clusterInfo\n");
	212	delete clusterInfo;
	213	printf("Closing esdFile\n");
	214	esdFile->Close();
	215	printf("Deleting esd\n");
	216	delete esd;
	217	// delete padInfo;
	218	cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
	219	}
	220
	221	//-----------------------------------------------------------------------
	222	void Prepare()
	223	{
	224	/// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
	225
8468d8c6	226	gRandom->SetSeed(0);
	227
	228	// Import TGeo geometry (needed for track extrapolation)
fcefbac4	229	if (!gGeoManager) {
	230	AliGeomManager::LoadGeometry("geometry.root");
	231	if (!gGeoManager) {
	232	Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
	233	return;
	234	}
	235	}
	236
	237	// set mag field
	238	printf("Loading field map...\n");
	239	AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
	240	AliTracker::SetFieldMap(field, kFALSE);
fcefbac4	241
	242	// Load mapping
	243	AliCDBManager* man = AliCDBManager::Instance();
	244	man->SetDefaultStorage("local://$ALICE_ROOT");
	245	man->SetRun(0);
	246	if ( ! AliMpCDB::LoadDDLStore() ) {
	247	Error("MUONRefit","Could not access mapping from OCDB !");
	248	exit(-1);
	249	}
	250
8468d8c6	251	// eventually set reconstruction parameters for refit (otherwise read from OCDB)
8468d8c6	252	/* AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
fcefbac4	253	muonRecoParam->Print("FULL");
fcefbac4	254	AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
8468d8c6	255	*/
fcefbac4	256	}
	257
	258	//-----------------------------------------------------------------------
	259	TTree* GetESDTree(TFile *esdFile)
	260	{
	261	/// Check that the file is properly open
	262	/// Return pointer to the ESD Tree
	263
	264	if (!esdFile \|\| !esdFile->IsOpen()) {
	265	Error("GetESDTree", "opening ESD file failed");
	266	exit(-1);
	267	}
	268
	269	TTree* tree = (TTree*) esdFile->Get("esdTree");
	270	if (!tree) {
	271	Error("GetESDTree", "no ESD tree found");
	272	exit(-1);
	273	}
	274
	275	return tree;
	276
	277	}
	278