[u/mrichter/AliRoot.git] / MUON / AliMUONReconstructor.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 for MUON reconstruction                                             //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////
#include <TParticle.h>
#include "TArrayF.h"
#include "AliRunLoader.h"
#include "AliRun.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliESD.h"
#include "AliMUONData.h"
#include "AliMUONEventReconstructor.h"
#include "AliMUONClusterReconstructor.h"
#include "AliMUONTriggerDecision.h"
#include "AliMUONClusterFinderVS.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTriggerTrack.h"
#include "AliESDMuonTrack.h"
#include "AliMUONReconstructor.h"

ClassImp(AliMUONReconstructor)
//_____________________________________________________________________________
AliMUONReconstructor::AliMUONReconstructor()
{
}
//_____________________________________________________________________________
AliMUONReconstructor::~AliMUONReconstructor()
{
}
//_____________________________________________________________________________
void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
{
//  AliLoader
  AliLoader* loader = runLoader->GetLoader("MUONLoader");
  Int_t nEvents = runLoader->GetNumberOfEvents();

// used local container for each method
// passing fLoader as argument, could be avoided ???
  AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(loader);
  AliMUONData* dataEvent = recoEvent->GetMUONData();

  AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
  AliMUONData* dataCluster = recoCluster->GetMUONData();

  AliMUONTriggerDecision* trigDec = new AliMUONTriggerDecision(loader);
  AliMUONData* dataTrig = trigDec->GetMUONData();


  for (Int_t i = 0; i < 10; i++) {
    AliMUONClusterFinderVS *RecModel = new AliMUONClusterFinderVS();
    RecModel->SetGhostChi2Cut(10);
    recoCluster->SetReconstructionModel(i,RecModel);
  } 

  loader->LoadDigits("READ");
  loader->LoadRecPoints("RECREATE");
  loader->LoadTracks("RECREATE");
  
  //   Loop over events              
  for(Int_t ievent = 0; ievent < nEvents; ievent++) {
    printf("Event %d\n",ievent);
    runLoader->GetEvent(ievent);

    //----------------------- digit2cluster & Digits2Trigger -------------------
    if (!loader->TreeR()) loader->MakeRecPointsContainer();
     
    // tracking branch
    dataCluster->MakeBranch("RC");
    dataCluster->SetTreeAddress("D,RC");
    recoCluster->Digits2Clusters(); 
    dataCluster->Fill("RC"); 

    // trigger branch
    dataTrig->MakeBranch("GLT");
    dataTrig->SetTreeAddress("D,GLT");
    trigDec->Digits2Trigger(); 
    dataTrig->Fill("GLT");

    loader->WriteRecPoints("OVERWRITE");

    //---------------------------- Track & TriggerTrack ---------------------
    if (!loader->TreeT()) loader->MakeTracksContainer();

    // trigger branch
    dataEvent->MakeBranch("RL"); //trigger track
    dataEvent->SetTreeAddress("RL");
    recoEvent->EventReconstructTrigger();
    dataEvent->Fill("RL");

    // tracking branch
    dataEvent->MakeBranch("RT"); //track
    dataEvent->SetTreeAddress("RT");
    recoEvent->EventReconstruct();
    dataEvent->Fill("RT");

    loader->WriteTracks("OVERWRITE");  
  
    //--------------------------- Resetting branches -----------------------
    dataCluster->ResetDigits();
    dataCluster->ResetRawClusters();

    dataTrig->ResetDigits();
    dataTrig->ResetTrigger();

    dataEvent->ResetRawClusters();
    dataEvent->ResetTrigger();
    dataEvent->ResetRecTracks();
    dataEvent->ResetRecTriggerTracks();
  
  }
  loader->UnloadDigits();
  loader->UnloadRecPoints();
  loader->UnloadTracks();

  delete recoCluster;
  delete recoEvent;
  delete trigDec;
}
//_____________________________________________________________________________
void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
{
  TClonesArray* recTracksArray = 0;
  TClonesArray* recTrigTracksArray = 0;
  
  AliLoader* loader = runLoader->GetLoader("MUONLoader");
  loader->LoadTracks("READ");
  AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");

   // declaration  
  Int_t iEvent, nPart;
  Int_t nTrackHits, nPrimary;
  Double_t fitFmin;
  TArrayF vertex(3);

  Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
  Double_t xRec, yRec, zRec, chi2MatchTrigger;
  Bool_t matchTrigger;

  // setting pointer for tracks, triggertracks & trackparam at vertex
  AliMUONTrack* recTrack = 0;
  AliMUONTrackParam* trackParam = 0;
  AliMUONTriggerTrack* recTriggerTrack = 0;
  TParticle* particle = new TParticle();
  AliGenEventHeader* header = 0;
  iEvent = runLoader->GetEventNumber(); 
  runLoader->GetEvent(iEvent);

  // vertex calculation (maybe it exists already somewhere else)
  vertex[0] = vertex[1] = vertex[2] = 0.;
  nPrimary = 0;
  if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
    header->PrimaryVertex(vertex);
  } else {
    runLoader->LoadKinematics("READ");
    runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
    nPart = (Int_t)runLoader->TreeK()->GetEntries();
    for(Int_t iPart = 0; iPart < nPart; iPart++) {
      runLoader->TreeK()->GetEvent(iPart);
      if (particle->GetFirstMother() == -1) {
	vertex[0] += particle->Vx();
	vertex[1] += particle->Vy();
	vertex[2] += particle->Vz();
	nPrimary++;
      }
      if (nPrimary) {
	vertex[0] /= (double)nPrimary;
	vertex[1] /= (double)nPrimary;
	vertex[2] /= (double)nPrimary;
      }
    }
  }
  // setting ESD MUON class
  AliESDMuonTrack* ESDTrack = new  AliESDMuonTrack() ;

  //-------------------- trigger tracks-------------
  Long_t trigPat = 0;
  muonData->SetTreeAddress("RL");
  muonData->GetRecTriggerTracks();
  recTrigTracksArray = muonData->RecTriggerTracks();

  // ready global trigger pattern from first track
  if (recTrigTracksArray) 
    recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
  if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();

  //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
 
  // -------------------- tracks-------------
  muonData->SetTreeAddress("RT");
  muonData->GetRecTracks();
  recTracksArray = muonData->RecTracks();
        
  Int_t nRecTracks = 0;
  if (recTracksArray)
    nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
  
  // loop over tracks
  for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {

    // reading info from tracks
    recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);

    trackParam = recTrack->GetTrackParamAtVertex();

    bendingSlope            = trackParam->GetBendingSlope();
    nonBendingSlope         = trackParam->GetNonBendingSlope();
    inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
    xRec  = trackParam->GetNonBendingCoor();
    yRec  = trackParam->GetBendingCoor();
    zRec  = trackParam->GetZ();

    nTrackHits       = recTrack->GetNTrackHits();
    fitFmin          = recTrack->GetFitFMin();
    matchTrigger     = recTrack->GetMatchTrigger();
    chi2MatchTrigger = recTrack->GetChi2MatchTrigger();

    // setting data member of ESD MUON
    ESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
    ESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
    ESDTrack->SetThetaY(TMath::ATan(bendingSlope));
    ESDTrack->SetZ(vertex[2]);
    ESDTrack->SetBendingCoor(vertex[1]); // calculate vertex at ESD or Tracking level ?
    ESDTrack->SetNonBendingCoor(vertex[0]);
    ESDTrack->SetChi2(fitFmin);
    ESDTrack->SetNHit(nTrackHits);
    ESDTrack->SetMatchTrigger(matchTrigger);
    ESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);

    // storing ESD MUON Track into ESD Event 
    if (nRecTracks != 0)  
      esd->AddMuonTrack(ESDTrack);
  } // end loop tracks

  // add global trigger pattern
  if (nRecTracks != 0)  
    esd->SetTrigger(trigPat);

  // reset muondata
  muonData->ResetRecTracks();
  muonData->ResetRecTriggerTracks();

  //} // end loop on event  
  loader->UnloadTracks(); 
  if (!header)
    runLoader->UnloadKinematics();
  delete ESDTrack;
  delete muonData;
  // delete particle;
}
Commit	Line	Data
cac2eb58	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	///////////////////////////////////////////////////////////////////////////////
	19	// //
b002ac70	20	// class for MUON reconstruction //
cac2eb58	21	// //
cac2eb58	22	///////////////////////////////////////////////////////////////////////////////
b002ac70	23	#include <TParticle.h>
b002ac70	24	#include "TArrayF.h"
cac2eb58	25	#include "AliRunLoader.h"
cac2eb58	26	#include "AliRun.h"
b002ac70	27	#include "AliHeader.h"
b002ac70	28	#include "AliGenEventHeader.h"
cac2eb58	29	#include "AliESD.h"
	30	#include "AliMUONData.h"
	31	#include "AliMUONEventReconstructor.h"
	32	#include "AliMUONClusterReconstructor.h"
	33	#include "AliMUONTriggerDecision.h"
	34	#include "AliMUONClusterFinderVS.h"
	35	#include "AliMUONTrack.h"
	36	#include "AliMUONTrackParam.h"
	37	#include "AliMUONTriggerTrack.h"
	38	#include "AliESDMuonTrack.h"
	39	#include "AliMUONReconstructor.h"
	40
	41	ClassImp(AliMUONReconstructor)
	42	//_____________________________________________________________________________
	43	AliMUONReconstructor::AliMUONReconstructor()
8789635b	44	{
	45	}
	46	//_____________________________________________________________________________
	47	AliMUONReconstructor::~AliMUONReconstructor()
cac2eb58	48	{
	49	}
	50	//_____________________________________________________________________________
	51	void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
	52	{
	53	// AliLoader
	54	AliLoader* loader = runLoader->GetLoader("MUONLoader");
	55	Int_t nEvents = runLoader->GetNumberOfEvents();
	56
	57	// used local container for each method
	58	// passing fLoader as argument, could be avoided ???
	59	AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(loader);
	60	AliMUONData* dataEvent = recoEvent->GetMUONData();
	61
	62	AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
	63	AliMUONData* dataCluster = recoCluster->GetMUONData();
	64
	65	AliMUONTriggerDecision* trigDec = new AliMUONTriggerDecision(loader);
	66	AliMUONData* dataTrig = trigDec->GetMUONData();
	67
	68
	69	for (Int_t i = 0; i < 10; i++) {
	70	AliMUONClusterFinderVS *RecModel = new AliMUONClusterFinderVS();
	71	RecModel->SetGhostChi2Cut(10);
	72	recoCluster->SetReconstructionModel(i,RecModel);
	73	}
	74
	75	loader->LoadDigits("READ");
	76	loader->LoadRecPoints("RECREATE");
	77	loader->LoadTracks("RECREATE");
	78
	79	// Loop over events
	80	for(Int_t ievent = 0; ievent < nEvents; ievent++) {
	81	printf("Event %d\n",ievent);
	82	runLoader->GetEvent(ievent);
	83
	84	//----------------------- digit2cluster & Digits2Trigger -------------------
	85	if (!loader->TreeR()) loader->MakeRecPointsContainer();
	86
	87	// tracking branch
	88	dataCluster->MakeBranch("RC");
	89	dataCluster->SetTreeAddress("D,RC");
	90	recoCluster->Digits2Clusters();
	91	dataCluster->Fill("RC");
	92
	93	// trigger branch
	94	dataTrig->MakeBranch("GLT");
	95	dataTrig->SetTreeAddress("D,GLT");
	96	trigDec->Digits2Trigger();
	97	dataTrig->Fill("GLT");
	98
	99	loader->WriteRecPoints("OVERWRITE");
	100
	101	//---------------------------- Track & TriggerTrack ---------------------
	102	if (!loader->TreeT()) loader->MakeTracksContainer();
	103
	104	// trigger branch
	105	dataEvent->MakeBranch("RL"); //trigger track
	106	dataEvent->SetTreeAddress("RL");
	107	recoEvent->EventReconstructTrigger();
	108	dataEvent->Fill("RL");
	109
	110	// tracking branch
	111	dataEvent->MakeBranch("RT"); //track
112	dataEvent->SetTreeAddress("RT");
113	recoEvent->EventReconstruct();
114	dataEvent->Fill("RT");
115
116	loader->WriteTracks("OVERWRITE");
117
118	//--------------------------- Resetting branches -----------------------
119	dataCluster->ResetDigits();
120	dataCluster->ResetRawClusters();
121
122	dataTrig->ResetDigits();
123	dataTrig->ResetTrigger();
124
125	dataEvent->ResetRawClusters();
126	dataEvent->ResetTrigger();
127	dataEvent->ResetRecTracks();
128	dataEvent->ResetRecTriggerTracks();
129
130	}
131	loader->UnloadDigits();
132	loader->UnloadRecPoints();
133	loader->UnloadTracks();
134
135	delete recoCluster;
136	delete recoEvent;
137	delete trigDec;
138	}
139	//_____________________________________________________________________________
140	void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
141	{
142	TClonesArray* recTracksArray = 0;
143	TClonesArray* recTrigTracksArray = 0;
144
145	AliLoader* loader = runLoader->GetLoader("MUONLoader");
146	loader->LoadTracks("READ");
cac2eb58	147	AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
	148
	149	// declaration
b002ac70	150	Int_t iEvent, nPart;
b002ac70	151	Int_t nTrackHits, nPrimary;
cac2eb58	152	Double_t fitFmin;
b002ac70	153	TArrayF vertex(3);
cac2eb58	154
	155	Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
	156	Double_t xRec, yRec, zRec, chi2MatchTrigger;
	157	Bool_t matchTrigger;
	158
b002ac70	159	// setting pointer for tracks, triggertracks & trackparam at vertex
cac2eb58	160	AliMUONTrack* recTrack = 0;
	161	AliMUONTrackParam* trackParam = 0;
	162	AliMUONTriggerTrack* recTriggerTrack = 0;
b002ac70	163	TParticle* particle = new TParticle();
b002ac70	164	AliGenEventHeader* header = 0;
cac2eb58	165	iEvent = runLoader->GetEventNumber();
	166	runLoader->GetEvent(iEvent);
	167
b002ac70	168	// vertex calculation (maybe it exists already somewhere else)
	169	vertex[0] = vertex[1] = vertex[2] = 0.;
	170	nPrimary = 0;
	171	if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
	172	header->PrimaryVertex(vertex);
	173	} else {
	174	runLoader->LoadKinematics("READ");
	175	runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
	176	nPart = (Int_t)runLoader->TreeK()->GetEntries();
	177	for(Int_t iPart = 0; iPart < nPart; iPart++) {
	178	runLoader->TreeK()->GetEvent(iPart);
	179	if (particle->GetFirstMother() == -1) {
	180	vertex[0] += particle->Vx();
	181	vertex[1] += particle->Vy();
	182	vertex[2] += particle->Vz();
	183	nPrimary++;
	184	}
	185	if (nPrimary) {
	186	vertex[0] /= (double)nPrimary;
	187	vertex[1] /= (double)nPrimary;
	188	vertex[2] /= (double)nPrimary;
	189	}
	190	}
	191	}
cac2eb58	192	// setting ESD MUON class
	193	AliESDMuonTrack* ESDTrack = new AliESDMuonTrack() ;
	194
	195	//-------------------- trigger tracks-------------
	196	Long_t trigPat = 0;
	197	muonData->SetTreeAddress("RL");
	198	muonData->GetRecTriggerTracks();
	199	recTrigTracksArray = muonData->RecTriggerTracks();
	200
	201	// ready global trigger pattern from first track
	202	if (recTrigTracksArray)
	203	recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
	204	if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
	205
	206	//printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
	207
	208	// -------------------- tracks-------------
	209	muonData->SetTreeAddress("RT");
	210	muonData->GetRecTracks();
	211	recTracksArray = muonData->RecTracks();
	212
	213	Int_t nRecTracks = 0;
	214	if (recTracksArray)
	215	nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
	216
	217	// loop over tracks
	218	for (Int_t iRecTracks = 0; iRecTracks < nRecTracks; iRecTracks++) {
	219
	220	// reading info from tracks
	221	recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
	222
	223	trackParam = recTrack->GetTrackParamAtVertex();
	224
	225	bendingSlope = trackParam->GetBendingSlope();
	226	nonBendingSlope = trackParam->GetNonBendingSlope();
	227	inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
	228	xRec = trackParam->GetNonBendingCoor();
	229	yRec = trackParam->GetBendingCoor();
	230	zRec = trackParam->GetZ();
	231
	232	nTrackHits = recTrack->GetNTrackHits();
	233	fitFmin = recTrack->GetFitFMin();
	234	matchTrigger = recTrack->GetMatchTrigger();
	235	chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
	236
	237	// setting data member of ESD MUON
	238	ESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
	239	ESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
	240	ESDTrack->SetThetaY(TMath::ATan(bendingSlope));
b002ac70	241	ESDTrack->SetZ(vertex[2]);
	242	ESDTrack->SetBendingCoor(vertex[1]); // calculate vertex at ESD or Tracking level ?
	243	ESDTrack->SetNonBendingCoor(vertex[0]);
cac2eb58	244	ESDTrack->SetChi2(fitFmin);
	245	ESDTrack->SetNHit(nTrackHits);
	246	ESDTrack->SetMatchTrigger(matchTrigger);
	247	ESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
	248
	249	// storing ESD MUON Track into ESD Event
	250	if (nRecTracks != 0)
	251	esd->AddMuonTrack(ESDTrack);
	252	} // end loop tracks
	253
	254	// add global trigger pattern
	255	if (nRecTracks != 0)
	256	esd->SetTrigger(trigPat);
	257
	258	// reset muondata
	259	muonData->ResetRecTracks();
	260	muonData->ResetRecTriggerTracks();
	261
	262	//} // end loop on event
	263	loader->UnloadTracks();
b002ac70	264	if (!header)
b002ac70	265	runLoader->UnloadKinematics();
cac2eb58	266	delete ESDTrack;
cac2eb58	267	delete muonData;
b002ac70	268	// delete particle;
cac2eb58	269	}