[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 <TRandom.h>
#include <TROOT.h>
#include <TMath.h>

// STEER includes
#include "AliMagF.h"
#include "AliTracker.h"
#include "AliESDEvent.h"
#include "AliRecoParam.h"
#include "AliCDBManager.h"
#include "AliRunLoader.h"
#include "AliLoader.h"

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

const Int_t printLevel = 1;

void Prepare();
TTree* GetESDTree(TFile *esdFile);
UInt_t buildClusterMap(AliMUONTrack &track);

//-----------------------------------------------------------------------
void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root",
		     const char* inFileName = "galice.root", const char* outFileName = "clusterInfo.root")
{
  /// 1) if (esdFileName != "")
  /// loop over ESD event and fill AliMUONClusterInfo object with cluster + corresponding track parameters;
  /// 2) if (inFileName != "")
  /// loop over RecPoints and fill AliMUONClusterInfo object with cluster not attached to a track;
  /// 3) write results in a new root file.
  ///
  /// ******************************************* WARNING ******************************************* ///
  /// track parameters at each cluster are recomputed by the interface using Kalman filter + Smoother ///
  /// (It can be changed by resetting the tracker in the interface with a new recoParam object)       ///
  /// and the magnetic field set in the function prepare()                                            ///
  /// ******************************************* WARNING ******************************************* ///
  
  Bool_t useESD = (strcmp(esdFileName,""));
  Bool_t useRecPoints = (strcmp(inFileName,""));
  if (!useESD && !useRecPoints) {
    Error("MUONClusterInfo","you must provide ESD and/or galice root file(s)");
    return;
  }
  
  AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
  AliMUONPadInfo padInfo;
  AliMUONCalibrationData* calibData = 0x0;
  AliMUONESDInterface esdInterface;
  AliMUONVClusterStore* clusterStore = 0x0;
  AliMUONVDigitStore* digitStore = 0x0;
  
  // prepare the refitting during ESD->MUON conversion
  Prepare();
  
  // open the ESD file and tree and connect the ESD event
  TFile* esdFile = 0x0;
  TTree* esdTree = 0x0;
  AliESDEvent* esd = 0x0;
  if (useESD) {
    esdFile = TFile::Open(esdFileName);
    esdTree = GetESDTree(esdFile);
    esd = new AliESDEvent();
    esd->ReadFromTree(esdTree);
  }
  
  // get the cluster from RecPoints
  AliRunLoader * rl = 0x0;
  AliLoader* MUONLoader = 0x0;
  if (useRecPoints) {
    rl = AliRunLoader::Open(inFileName,"MUONLoader");
    MUONLoader = rl->GetDetectorLoader("MUON");
    MUONLoader->LoadRecPoints("READ");   
    MUONLoader->LoadDigits("READ");
  }
  
  // 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 events
  if (useESD) {
    if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
    else nevents = (Int_t)esdTree->GetEntries();
  } else {
    if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)rl->GetNumberOfEvents());
    else nevents = (Int_t)rl->GetNumberOfEvents();
  }
  for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
    
    //----------------------------------------------//
    // -------------- process event --------------- //
    //----------------------------------------------//
    // get the ESD of current event
    if (useESD) {
      esdTree->GetEvent(iEvent);
      if (!esd) {
        Error("MUONClusterInfo", "no ESD object found for event %d", iEvent);
        return;
      }
      // load the current esd event
      esdInterface.LoadEvent(*esd);
      // get digit store
      if (!useRecPoints) digitStore = esdInterface.GetDigits();
    }
    
    if (useRecPoints) {
      if (!(rl->GetEvent(iEvent) == 0)) {
        Error("MUONClusterInfo", "unable to load event %d", iEvent);
	return;
      }
      // get the clusters of current event
      TTree* treeR = MUONLoader->TreeR();
      clusterStore = AliMUONVClusterStore::Create(*treeR);
      if ( clusterStore != 0x0 ) {
	clusterStore->Clear();
	clusterStore->Connect(*treeR);
	treeR->GetEvent(0);
      }
      // get the digits of current event
      TTree* treeD = MUONLoader->TreeD();
      digitStore = AliMUONVDigitStore::Create(*treeD);
      if ( digitStore != 0x0 ) {
	digitStore->Clear();
	digitStore->Connect(*treeD);
	treeD->GetEvent(0);
      }
    }
    
    // prepare access to calibration data
    if (useESD && !calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
    else if (!calibData) calibData = new AliMUONCalibrationData(rl->GetRunNumber());
    
    //----------------------------------------------//
    // ------------- fill cluster info ------------ //
    //----------------------------------------------//
    // --- loop over the refitted tracks ---
    if (useESD) {
      
      TIter nextTrack(esdInterface.CreateTrackIterator());
      AliMUONTrack* track;
      while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
	
	UInt_t muonClusterMap = buildClusterMap(*track);
	
	// loop over clusters
	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
	while (trackParam) {
	  clusterInfo->Clear("C");
	  
	  // fill cluster info
	  AliMUONVCluster* cluster = trackParam->GetClusterPtr();
	  clusterInfo->SetRunId(esd->GetRunNumber());
	  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->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
	  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());
	  clusterInfo->SetTrackNHits(track->GetNClusters());
	  clusterInfo->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 = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
	    AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
	    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());
	    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);
	  }
	  
	  // remove clusters attached to a track
	  if (useRecPoints) {
	    AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID());
	    if (cl) clusterStore->Remove(*cl);
	  }
	  
	  // fill cluster info tree
	  clusterInfoTree->Fill();
	  
	  trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
	}
	
      }
      
    }
    
    // --- loop over clusters not attached to a track ---
    if (useRecPoints) {
      
      TIter nextCluster(clusterStore->CreateIterator());
      AliMUONVCluster *cluster;
      while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) {
	
	clusterInfo->Clear("C");
	
	// fill cluster info
	clusterInfo->SetRunId(rl->GetRunNumber());
	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 dummy track info
	clusterInfo->SetTrackId(0);
	clusterInfo->SetTrackXY(0.,0.);
	clusterInfo->SetTrackThetaXY(0.,0.);
	clusterInfo->SetTrackP(0.);
	clusterInfo->SetTrackXYErr(0.,0.);
	clusterInfo->SetTrackChi2(0.);
	clusterInfo->SetTrackCharge(0);
	clusterInfo->SetTrackNHits(0);
	clusterInfo->SetTrackChamberHitMap(0);
	
	// 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 = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
	  AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
	  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());
	  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();
	
      }
      
      delete digitStore;
      delete clusterStore;
      
    }
    
  }
  
  // ...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;
  if (useRecPoints) {
    MUONLoader->UnloadDigits();
    MUONLoader->UnloadRecPoints();
    delete rl;
  }
  if (useESD) {
    esdFile->Close();
    delete esd;
  }
  cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
}

//-----------------------------------------------------------------------
void Prepare()
{
  /// Set the magnetic field, the mapping and the reconstruction parameters
  
  gRandom->SetSeed(0);
  
  // set mag field
  if (!TGeoGlobalMagField::Instance()->GetField()) {
    printf("Loading field map...\n");
    AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
    TGeoGlobalMagField::Instance()->SetField(field);
  }
  
  // Load mapping
  AliCDBManager* man = AliCDBManager::Instance();
  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  man->SetRun(0);
  if ( ! AliMpCDB::LoadDDLStore() ) {
    Error("Prepare","Could not access mapping from OCDB !");
    exit(-1);
  }
  
  // Reset the reconstruction parameters for track refitting if needed
  // (by default will use Kalman filter + Smoother)
  //  AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
  //  AliMUONESDInterface::ResetTracker(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;
  
}

      
//-----------------------------------------------------------------------
UInt_t 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
fcefbac4	1	/**************************************************************************
41c52850	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	**************************************************************************/
fcefbac4	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>
fcefbac4	31	#include <TRandom.h>
fcefbac4	32	#include <TROOT.h>
41c52850	33	#include <TMath.h>
fcefbac4	34
fcefbac4	35	// STEER includes
f7a1cc68	36	#include "AliMagF.h"
fcefbac4	37	#include "AliTracker.h"
fcefbac4	38	#include "AliESDEvent.h"
fcefbac4	39	#include "AliRecoParam.h"
fcefbac4	40	#include "AliCDBManager.h"
41c52850	41	#include "AliRunLoader.h"
41c52850	42	#include "AliLoader.h"
fcefbac4	43
fcefbac4	44	// MUON includes
41c52850	45	#include "AliMpConstants.h"
fcefbac4	46	#include "AliMpCDB.h"
	47	#include "AliMpSegmentation.h"
	48	#include "AliMpVSegmentation.h"
	49	#include "AliMpPad.h"
	50	#include "AliMUONCalibrationData.h"
	51	#include "AliMUONVCalibParam.h"
	52	#include "AliMUONPadInfo.h"
	53	#include "AliMUONClusterInfo.h"
	54	#include "AliMUONRecoParam.h"
	55	#include "AliMUONESDInterface.h"
fcefbac4	56	#include "AliMUONVDigit.h"
	57	#include "AliMUONVDigitStore.h"
	58	#include "AliMUONVCluster.h"
41c52850	59	#include "AliMUONVClusterStore.h"
fcefbac4	60	#include "AliMUONTrack.h"
fcefbac4	61	#include "AliMUONTrackParam.h"
fcefbac4	62	#endif
	63
	64	const Int_t printLevel = 1;
	65
	66	void Prepare();
	67	TTree* GetESDTree(TFile *esdFile);
41c52850	68	UInt_t buildClusterMap(AliMUONTrack &track);
fcefbac4	69
fcefbac4	70	//-----------------------------------------------------------------------
41c52850	71	void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root",
41c52850	72	const char* inFileName = "galice.root", const char* outFileName = "clusterInfo.root")
fcefbac4	73	{
41c52850	74	/// 1) if (esdFileName != "")
	75	/// loop over ESD event and fill AliMUONClusterInfo object with cluster + corresponding track parameters;
	76	/// 2) if (inFileName != "")
	77	/// loop over RecPoints and fill AliMUONClusterInfo object with cluster not attached to a track;
	78	/// 3) write results in a new root file.
	79	///
	80	/// ***************************************** WARNING ***************************************** ///
	81	/// track parameters at each cluster are recomputed by the interface using Kalman filter + Smoother ///
	82	/// (It can be changed by resetting the tracker in the interface with a new recoParam object) ///
	83	/// and the magnetic field set in the function prepare() ///
	84	/// ***************************************** WARNING ***************************************** ///
	85
	86	Bool_t useESD = (strcmp(esdFileName,""));
	87	Bool_t useRecPoints = (strcmp(inFileName,""));
	88	if (!useESD && !useRecPoints) {
	89	Error("MUONClusterInfo","you must provide ESD and/or galice root file(s)");
	90	return;
	91	}
fcefbac4	92
	93	AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
	94	AliMUONPadInfo padInfo;
	95	AliMUONCalibrationData* calibData = 0x0;
8468d8c6	96	AliMUONESDInterface esdInterface;
41c52850	97	AliMUONVClusterStore* clusterStore = 0x0;
41c52850	98	AliMUONVDigitStore* digitStore = 0x0;
fcefbac4	99
8468d8c6	100	// prepare the refitting during ESD->MUON conversion
fcefbac4	101	Prepare();
fcefbac4	102
fcefbac4	103	// open the ESD file and tree and connect the ESD event
41c52850	104	TFile* esdFile = 0x0;
	105	TTree* esdTree = 0x0;
	106	AliESDEvent* esd = 0x0;
	107	if (useESD) {
	108	esdFile = TFile::Open(esdFileName);
	109	esdTree = GetESDTree(esdFile);
	110	esd = new AliESDEvent();
	111	esd->ReadFromTree(esdTree);
	112	}
	113
	114	// get the cluster from RecPoints
	115	AliRunLoader * rl = 0x0;
	116	AliLoader* MUONLoader = 0x0;
	117	if (useRecPoints) {
	118	rl = AliRunLoader::Open(inFileName,"MUONLoader");
	119	MUONLoader = rl->GetDetectorLoader("MUON");
	120	MUONLoader->LoadRecPoints("READ");
	121	MUONLoader->LoadDigits("READ");
	122	}
fcefbac4	123
	124	// prepare the output tree
	125	gROOT->cd();
	126	TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
	127	clusterInfoFile->SetCompressionLevel(1);
	128
	129	TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
	130	clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
	131
	132	// timer start...
	133	TStopwatch timer;
	134
41c52850	135	// Loop over events
	136	if (useESD) {
	137	if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
	138	else nevents = (Int_t)esdTree->GetEntries();
	139	} else {
	140	if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)rl->GetNumberOfEvents());
	141	else nevents = (Int_t)rl->GetNumberOfEvents();
	142	}
fcefbac4	143	for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
	144
	145	//----------------------------------------------//
	146	// -------------- process event --------------- //
	147	//----------------------------------------------//
	148	// get the ESD of current event
41c52850	149	if (useESD) {
	150	esdTree->GetEvent(iEvent);
	151	if (!esd) {
	152	Error("MUONClusterInfo", "no ESD object found for event %d", iEvent);
	153	return;
	154	}
	155	// load the current esd event
	156	esdInterface.LoadEvent(*esd);
	157	// get digit store
	158	if (!useRecPoints) digitStore = esdInterface.GetDigits();
fcefbac4	159	}
fcefbac4	160
41c52850	161	if (useRecPoints) {
	162	if (!(rl->GetEvent(iEvent) == 0)) {
	163	Error("MUONClusterInfo", "unable to load event %d", iEvent);
	164	return;
	165	}
	166	// get the clusters of current event
	167	TTree* treeR = MUONLoader->TreeR();
	168	clusterStore = AliMUONVClusterStore::Create(*treeR);
	169	if ( clusterStore != 0x0 ) {
	170	clusterStore->Clear();
	171	clusterStore->Connect(*treeR);
	172	treeR->GetEvent(0);
	173	}
	174	// get the digits of current event
	175	TTree* treeD = MUONLoader->TreeD();
	176	digitStore = AliMUONVDigitStore::Create(*treeD);
	177	if ( digitStore != 0x0 ) {
	178	digitStore->Clear();
	179	digitStore->Connect(*treeD);
	180	treeD->GetEvent(0);
	181	}
	182	}
fcefbac4	183
41c52850	184	// prepare access to calibration data
	185	if (useESD && !calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
	186	else if (!calibData) calibData = new AliMUONCalibrationData(rl->GetRunNumber());
fcefbac4	187
fcefbac4	188	//----------------------------------------------//
	189	// ------------- fill cluster info ------------ //
	190	//----------------------------------------------//
41c52850	191	// --- loop over the refitted tracks ---
	192	if (useESD) {
	193
	194	TIter nextTrack(esdInterface.CreateTrackIterator());
	195	AliMUONTrack* track;
	196	while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
	197
	198	UInt_t muonClusterMap = buildClusterMap(*track);
	199
	200	// loop over clusters
	201	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
	202	while (trackParam) {
	203	clusterInfo->Clear("C");
	204
	205	// fill cluster info
	206	AliMUONVCluster* cluster = trackParam->GetClusterPtr();
	207	clusterInfo->SetRunId(esd->GetRunNumber());
	208	clusterInfo->SetEventId(iEvent);
	209	clusterInfo->SetZ(cluster->GetZ());
	210	clusterInfo->SetClusterId(cluster->GetUniqueID());
	211	clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
	212	clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
	213	clusterInfo->SetClusterChi2(cluster->GetChi2());
	214	clusterInfo->SetClusterCharge(cluster->GetCharge());
	215
	216	// fill track info
	217	clusterInfo->SetTrackId(track->GetUniqueID());
	218	clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
	219	clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
	220	clusterInfo->SetTrackP(trackParam->P());
	221	const TMatrixD paramCov = trackParam->GetCovariances();
	222	clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
	223	clusterInfo->SetTrackChi2(track->GetNormalizedChi2());
	224	clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
	225	clusterInfo->SetTrackNHits(track->GetNClusters());
	226	clusterInfo->SetTrackChamberHitMap(muonClusterMap);
	227
	228	// fill pad info if available
	229	for (Int_t i=0; i<cluster->GetNDigits(); i++) {
	230	AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
	231	if (!digit) continue;
	232
	233	// pad location
	234	const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
	235	GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
	236	AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
	237
	238	// calibration parameters
	239	AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
	240	AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
	241	Int_t manuChannel = digit->ManuChannel();
	242	Int_t planeType = 0;
	243	if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
	244	planeType = 1;
	245	}
	246
	247	// fill pad info
	248	padInfo.SetPadId(digit->GetUniqueID());
	249	padInfo.SetPadPlaneType(planeType);
	250	padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
	251	padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
	252	padInfo.SetPadCharge((Double_t)digit->Charge());
	253	padInfo.SetPadADC(digit->ADC());
	254	padInfo.SetSaturated(digit->IsSaturated());
255	padInfo.SetCalibrated(digit->IsCalibrated());
256	padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
257	padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
258	gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
259
260	clusterInfo->AddPad(padInfo);
261	}
262
263	// remove clusters attached to a track
264	if (useRecPoints) {
265	AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID());
266	if (cl) clusterStore->Remove(*cl);
267	}
268
269	// fill cluster info tree
270	clusterInfoTree->Fill();
271
272	trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
273	}
274
275	}
276
277	}
278
279	// --- loop over clusters not attached to a track ---
280	if (useRecPoints) {
fcefbac4	281
41c52850	282	TIter nextCluster(clusterStore->CreateIterator());
	283	AliMUONVCluster *cluster;
	284	while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) {
	285
fcefbac4	286	clusterInfo->Clear("C");
	287
	288	// fill cluster info
41c52850	289	clusterInfo->SetRunId(rl->GetRunNumber());
fcefbac4	290	clusterInfo->SetEventId(iEvent);
	291	clusterInfo->SetZ(cluster->GetZ());
	292	clusterInfo->SetClusterId(cluster->GetUniqueID());
	293	clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
	294	clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
	295	clusterInfo->SetClusterChi2(cluster->GetChi2());
	296	clusterInfo->SetClusterCharge(cluster->GetCharge());
	297
41c52850	298	// fill dummy track info
	299	clusterInfo->SetTrackId(0);
	300	clusterInfo->SetTrackXY(0.,0.);
	301	clusterInfo->SetTrackThetaXY(0.,0.);
	302	clusterInfo->SetTrackP(0.);
	303	clusterInfo->SetTrackXYErr(0.,0.);
	304	clusterInfo->SetTrackChi2(0.);
	305	clusterInfo->SetTrackCharge(0);
	306	clusterInfo->SetTrackNHits(0);
	307	clusterInfo->SetTrackChamberHitMap(0);
	308
fcefbac4	309	// fill pad info if available
	310	for (Int_t i=0; i<cluster->GetNDigits(); i++) {
	311	AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
	312	if (!digit) continue;
	313
	314	// pad location
	315	const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
41c52850	316	GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
168e9c4d	317	AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
fcefbac4	318
	319	// calibration parameters
	320	AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
	321	AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
	322	Int_t manuChannel = digit->ManuChannel();
41c52850	323	Int_t planeType = 0;
	324	if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
	325	planeType = 1;
	326	}
fcefbac4	327
	328	// fill pad info
	329	padInfo.SetPadId(digit->GetUniqueID());
41c52850	330	padInfo.SetPadPlaneType(planeType);
68828de2	331	padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
68828de2	332	padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
fcefbac4	333	padInfo.SetPadCharge((Double_t)digit->Charge());
	334	padInfo.SetPadADC(digit->ADC());
	335	padInfo.SetSaturated(digit->IsSaturated());
	336	padInfo.SetCalibrated(digit->IsCalibrated());
	337	padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
	338	padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
	339	gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
	340
	341	clusterInfo->AddPad(padInfo);
	342	}
	343
	344	// fill cluster info tree
	345	clusterInfoTree->Fill();
	346
fcefbac4	347	}
fcefbac4	348
41c52850	349	delete digitStore;
	350	delete clusterStore;
	351
fcefbac4	352	}
fcefbac4	353
fcefbac4	354	}
	355
	356	// ...timer stop
	357	timer.Stop();
	358	printf("Writing Tree\n");
	359	// write output tree
	360	clusterInfoFile->cd();
	361	clusterInfoTree->Write();
	362	printf("Deleting Tree\n");
	363	delete clusterInfoTree;
	364	printf("Closing File\n");
	365	clusterInfoFile->Close();
	366
	367	// free memory
	368	printf("Deleting calibData\n");
	369	delete calibData;
	370	printf("Deleting clusterInfo\n");
	371	delete clusterInfo;
41c52850	372	if (useRecPoints) {
	373	MUONLoader->UnloadDigits();
	374	MUONLoader->UnloadRecPoints();
	375	delete rl;
	376	}
	377	if (useESD) {
	378	esdFile->Close();
	379	delete esd;
	380	}
fcefbac4	381	cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
	382	}
	383
	384	//-----------------------------------------------------------------------
	385	void Prepare()
	386	{
41c52850	387	/// Set the magnetic field, the mapping and the reconstruction parameters
fcefbac4	388
8468d8c6	389	gRandom->SetSeed(0);
8468d8c6	390
fcefbac4	391	// set mag field
f7a1cc68	392	if (!TGeoGlobalMagField::Instance()->GetField()) {
	393	printf("Loading field map...\n");
	394	AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
	395	TGeoGlobalMagField::Instance()->SetField(field);
	396	}
fcefbac4	397
	398	// Load mapping
	399	AliCDBManager* man = AliCDBManager::Instance();
162637e4	400	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
fcefbac4	401	man->SetRun(0);
fcefbac4	402	if ( ! AliMpCDB::LoadDDLStore() ) {
41c52850	403	Error("Prepare","Could not access mapping from OCDB !");
fcefbac4	404	exit(-1);
	405	}
	406
ad3c6eda	407	// Reset the reconstruction parameters for track refitting if needed
ad3c6eda	408	// (by default will use Kalman filter + Smoother)
41c52850	409	// AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
41c52850	410	// AliMUONESDInterface::ResetTracker(muonRecoParam);
ad3c6eda	411
fcefbac4	412	}
	413
	414	//-----------------------------------------------------------------------
	415	TTree* GetESDTree(TFile *esdFile)
	416	{
	417	/// Check that the file is properly open
	418	/// Return pointer to the ESD Tree
	419
	420	if (!esdFile \|\| !esdFile->IsOpen()) {
	421	Error("GetESDTree", "opening ESD file failed");
	422	exit(-1);
	423	}
	424
	425	TTree* tree = (TTree*) esdFile->Get("esdTree");
	426	if (!tree) {
	427	Error("GetESDTree", "no ESD tree found");
	428	exit(-1);
	429	}
	430
	431	return tree;
	432
	433	}
	434
41c52850	435
	436	//-----------------------------------------------------------------------
	437	UInt_t buildClusterMap(AliMUONTrack &track)
	438	{
	439	/// Build the map of clusters in tracking chambers
	440
	441	UInt_t muonClusterMap = 0;
	442
	443	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
	444	while (trackParam) {
	445
	446	muonClusterMap \|= BIT(trackParam->GetClusterPtr()->GetChamberId());
	447
	448	trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
	449	}
	450
	451	return muonClusterMap;
	452
	453	}
	454