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

/**************************************************************************
 * Copyright(c) 1998-2000, 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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.9  2002/03/13 07:04:11  jchudoba
Connect only MUON branches when reading the event to speed up digitisation.

Revision 1.8  2002/02/22 12:14:21  morsch
Validate pad hit before digitization.

Revision 1.7  2001/11/22 11:15:41  jchudoba
Proper deletion of arrays (thanks to Rene Brun)

Revision 1.6  2001/11/02 12:55:45  jchudoba
cleanup of the code, add const to Get methods

Revision 1.5  2001/10/31 16:35:07  jchudoba
some functionality move to AliMUONTransientDigit class

Revision 1.4  2001/03/20 13:36:11  egangler
TFile memory leak and "too many files open" problem solved (?):
the backgroud file is open only once forever, and not once for each event.

Revision 1.3  2001/03/05 23:57:44  morsch
Writing of digit tree moved to macro.

Revision 1.2  2001/03/05 08:40:25  morsch
Method SortTracks(..) imported from AliMUON.

Revision 1.1  2001/02/02 14:11:53  morsch
AliMUONMerger prototype to be called by the merge manager.

*/

#include <TDirectory.h>
#include <TFile.h>
#include <TObjArray.h>
#include <TPDGCode.h>
#include <TTree.h> 

#include "AliHitMap.h"
#include "AliMUON.h"
#include "AliMUONChamber.h"
#include "AliMUONConstants.h"
#include "AliMUONDigit.h"
#include "AliMUONHit.h"
#include "AliMUONHitMapA1.h"
#include "AliMUONMerger.h"
#include "AliMUONPadHit.h"
#include "AliMUONTransientDigit.h"
#include "AliRun.h"

ClassImp(AliMUONMerger)

//___________________________________________
AliMUONMerger::AliMUONMerger()
{
// Default constructor    
    fEvNrSig = 0;
    fEvNrBgr = 0;
    fMerge   = kDigitize;
    fFnBgr   = 0;
    fHitMap  = 0;
    fList    = 0;
    fTrH1    = 0;
    fHitsBgr = 0;
    fPadHitsBgr = 0;
    fHitMap     = 0;
    fList       = 0;
    fBgrFile    = 0;
}

//------------------------------------------------------------------------
AliMUONMerger::~AliMUONMerger()
{
// Destructor
    if (fTrH1)       delete fTrH1;
    if (fHitsBgr)    delete fHitsBgr;
    if (fPadHitsBgr) delete fPadHitsBgr;
    if (fHitMap)     delete [] fHitMap;
    if (fList)       delete fList;
    if (fBgrFile)    delete fBgrFile;
}

//------------------------------------------------------------------------
Bool_t AliMUONMerger::Exists(const AliMUONPadHit *padhit) const
{
// test if the given padhit was already fired
    return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
}

//------------------------------------------------------------------------
void AliMUONMerger::Update(AliMUONPadHit *padhit)
{
// add new contribution to the fired padhit
    AliMUONTransientDigit *pdigit = 
      static_cast<AliMUONTransientDigit*>(
      fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));

    // update charge
    //
    Int_t iqpad    = padhit->QPad();        // charge per pad
    pdigit->AddSignal(iqpad);
    pdigit->AddPhysicsSignal(iqpad);		

    // update list of tracks
    //
    Int_t track, charge;    
    if (fSignal) {
      track = fTrack;
      charge = iqpad;
    } else {
      track = kBgTag;
      charge = kBgTag;
    }
    pdigit->UpdateTrackList(track,charge);
}

//------------------------------------------------------------------------
void AliMUONMerger::CreateNew(AliMUONPadHit *padhit)
{
// add new transient digit to the list, update hit map
    fList->AddAtAndExpand(
	new AliMUONTransientDigit(fNch,fDigits),fCounter);
    fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
    AliMUONTransientDigit* pdigit = 
      static_cast<AliMUONTransientDigit*>(fList->Last());
    // list of tracks
    Int_t track, charge;
    if (fSignal) {
      track = fTrack;
      charge = padhit->QPad();
    } else {
      track = kBgTag;
      charge = kBgTag;
    }
    pdigit->AddToTrackList(track,charge);
    fCounter++;
}


//------------------------------------------------------------------------
void AliMUONMerger::Init()
{
// Initialisation
    // open only once the background file !!
    if (fMerge && !fBgrFile) fBgrFile = InitBgr();
}


//------------------------------------------------------------------------
TFile* AliMUONMerger::InitBgr()
{
// Initialise background event
    TFile *file = new TFile(fFnBgr);
// add error checking later
    printf("\n AliMUONMerger has opened %s file with background event \n", fFnBgr);
    fHitsBgr     = new TClonesArray("AliMUONHit",1000);
    fPadHitsBgr  = new TClonesArray("AliMUONPadHit",1000);
    return file;
}

//------------------------------------------------------------------------
void AliMUONMerger::Digitise()
{

    // keep galice.root for signal and name differently the file for 
    // background when add! otherwise the track info for signal will be lost !
  
    AliMUONChamber*   iChamber;
    AliSegmentation*  segmentation;

    fList = new TObjArray;

    AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
    fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
    if (fMerge ) {
	fBgrFile->cd();
        //
	// Get Hits Tree header from file
	//if(fHitsBgr) fHitsBgr->Clear();     // Useless because line 327
	//if(fPadHitsBgr) fPadHitsBgr->Clear(); // Useless because line 328
	if(fTrH1) delete fTrH1;
	fTrH1 = 0;
	
	char treeName[20];
	sprintf(treeName,"TreeH%d",fEvNrBgr);
	fTrH1 = (TTree*)gDirectory->Get(treeName);
	if (!fTrH1) {
	    printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
	}
        //
	// Set branch addresses
	TBranch *branch;
	char branchname[20];
	sprintf(branchname,"%s",pMUON->GetName());
	if (fTrH1 && fHitsBgr) {
	    branch = fTrH1->GetBranch(branchname);
	    if (branch) branch->SetAddress(&fHitsBgr);
	}
	if (fTrH1 && fPadHitsBgr) {
	    branch = fTrH1->GetBranch("MUONCluster");
	    if (branch) branch->SetAddress(&fPadHitsBgr);
	}
    }
    //
    // loop over cathodes
    //
    fSignal = kTRUE;
    for (int icat = 0; icat < 2; icat++) { 
	fCounter = 0;
	for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
	    iChamber = &(pMUON->Chamber(i));
	    if (iChamber->Nsec() == 1 && icat == 1) {
		continue;
	    } else {
		segmentation = iChamber->SegmentationModel(icat+1);
	    }
	    fHitMap[i] = new AliMUONHitMapA1(segmentation, fList);
	}

//
//   Loop over tracks
//

	TTree *treeH  = gAlice->TreeH();
	Int_t ntracks = (Int_t) treeH->GetEntries();
	treeH->SetBranchStatus("*",0); // switch off all branches
        treeH->SetBranchStatus("MUON*",1); // switch on only MUON

	for (fTrack = 0; fTrack < ntracks; fTrack++) {
	    gAlice->ResetHits();
	    treeH->GetEntry(fTrack,0);
//
//   Loop over hits
	    for(AliMUONHit* mHit = (AliMUONHit*)pMUON->FirstHit(-1); 
		mHit;
		mHit = (AliMUONHit*)pMUON->NextHit()) 
	    {
		fNch = mHit->Chamber()-1;  // chamber number
		if (fNch > AliMUONConstants::NCh()-1) continue;
		iChamber = &(pMUON->Chamber(fNch));
		
//
// Loop over pad hits
		for (AliMUONPadHit* mPad =
			 (AliMUONPadHit*)pMUON->FirstPad(mHit,pMUON->PadHits());
		     mPad;
		     mPad = (AliMUONPadHit*)pMUON->NextPad(pMUON->PadHits()))
		{
		    Int_t cathode  = mPad->Cathode();      // cathode number
		    if (cathode != (icat+1)) continue;
		    Int_t iqpad    = Int_t(mPad->QPad());  // charge per pad
//		    segmentation = iChamber->SegmentationModel(cathode);
		    fDigits[0] = mPad->PadX();  
		    fDigits[1] = mPad->PadY();
		    if (!(fHitMap[fNch]->ValidateHit(fDigits[0], fDigits[1]))) continue;
		    fDigits[2] = icat;
		    fDigits[3] = iqpad;
		    fDigits[4] = iqpad;
		    if (mHit->Particle() == kMuonPlus ||
			mHit->Particle() == kMuonMinus) {
			fDigits[5] = mPad->HitNumber();
		    } else fDigits[5] = -1;

		    // build the list of fired pads and update the info

		    if (!Exists(mPad)) {
			CreateNew(mPad);
		    } else {
			Update(mPad);
		    } //  end if pdigit
		} //end loop over clusters
	    } // hit loop
	} // track loop

	// open the file with background
       
	if (fMerge) {
            fSignal = kFALSE;
	    ntracks = (Int_t)fTrH1->GetEntries();
//
//   Loop over tracks
//
	    for (fTrack = 0; fTrack < ntracks; fTrack++) {

		if (fHitsBgr)       fHitsBgr->Clear();
		if (fPadHitsBgr)    fPadHitsBgr->Clear();

		fTrH1->GetEvent(fTrack);
//
//   Loop over hits
		AliMUONHit* mHit;
		for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i) 
		{	
		    mHit   = (AliMUONHit*) (*fHitsBgr)[i];
		    fNch   = mHit->Chamber()-1;  // chamber number
		    iChamber = &(pMUON->Chamber(fNch));
//
// Loop over pad hits
		    for (AliMUONPadHit* mPad =
			     (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHitsBgr);
			 mPad;
			 mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHitsBgr))
		    {
			Int_t cathode  = mPad->Cathode();     // cathode number
			Int_t ipx      = mPad->PadX();        // pad number on X
			Int_t ipy      = mPad->PadY();        // pad number on Y
			Int_t iqpad    = Int_t(mPad->QPad()); // charge per pad
			if (!(fHitMap[fNch]->ValidateHit(ipx, ipy))) continue;

			if (cathode != (icat+1)) continue;
			fDigits[0] = ipx;
			fDigits[1] = ipy;
			fDigits[2] = icat;
			fDigits[3] = iqpad;
			fDigits[4] = 0;
			fDigits[5] = -1;
			
			// build the list of fired pads and update the info
			if (!Exists(mPad)) {
			    CreateNew(mPad);
			} else {
			    Update(mPad);
			} //  end if !Exists
		    } //end loop over clusters
		} // hit loop
	    } // track loop

	    TTree *treeK = gAlice->TreeK();
            TFile *file = NULL;
	    
	    if (treeK) file = treeK->GetCurrentFile();
	    file->cd();
	} // if fMerge

	Int_t tracks[kMAXTRACKS];
	Int_t charges[kMAXTRACKS];
	Int_t nentries = fList->GetEntriesFast();
	
	for (Int_t nent = 0; nent < nentries; nent++) {
	    AliMUONTransientDigit *address = (AliMUONTransientDigit*)fList->At(nent);
	    if (address == 0) continue; 
	    Int_t ich = address->Chamber();
	    Int_t   q = address->Signal(); 
	    iChamber = &(pMUON->Chamber(ich));
//
//  Digit Response (noise, threshold, saturation, ...)
	    AliMUONResponse * response = iChamber->ResponseModel();
	    q = response->DigitResponse(q);
	    
	    if (!q) continue;
	    
	    fDigits[0] = address->PadX();
	    fDigits[1] = address->PadY();
	    fDigits[2] = address->Cathode();
	    fDigits[3] = q;
	    fDigits[4] = address->Physics();
	    fDigits[5] = address->Hit();
	    
	    Int_t nptracks = address->GetNTracks();

	    if (nptracks > kMAXTRACKS) {
		printf("\n Attention - nptracks > kMAXTRACKS %d \n", nptracks);
		nptracks = kMAXTRACKS;
	    }
	    if (nptracks > 2) {
		printf("Attention - nptracks > 2  %d \n",nptracks);
		printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
	    }
	    for (Int_t tr = 0; tr < nptracks; tr++) {
		tracks[tr]   = address->GetTrack(tr);
		charges[tr]  = address->GetCharge(tr);
	    }      //end loop over list of tracks for one pad
            // Sort list of tracks according to charge
	    if (nptracks > 1) {
		SortTracks(tracks,charges,nptracks);
	    }
	    if (nptracks < kMAXTRACKS ) {
		for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
		    tracks[i]  = 0;
		    charges[i] = 0;
		}
	    }
	    
	    // fill digits
	    pMUON->AddDigits(ich,tracks,charges,fDigits);
	}
	gAlice->TreeD()->Fill();
	pMUON->ResetDigits();
	fList->Delete();

	
	for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
	    if (fHitMap[ii]) {
		delete fHitMap[ii];
		fHitMap[ii] = 0;
	    }
	}
    } //end loop over cathodes
    delete [] fHitMap;
    delete fList;
    
//  no need to delete ... and it makes a crash also
//    if (fHitsBgr)    fHitsBgr->Delete();
//    if (fPadHitsBgr) fPadHitsBgr->Delete();
    // gObjectTable->Print();
}


void AliMUONMerger::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
{
  //
  // Sort the list of tracks contributing to a given digit
  // Only the 3 most significant tracks are acctually sorted
  //
  
  //
  //  Loop over signals, only 3 times
  //
  
  Int_t qmax;
  Int_t jmax;
  Int_t idx[3] = {-2,-2,-2};
  Int_t jch[3] = {-2,-2,-2};
  Int_t jtr[3] = {-2,-2,-2};
  Int_t i,j,imax;
  
  if (ntr<3) imax=ntr;
  else imax=3;
  for(i=0;i<imax;i++){
    qmax=0;
    jmax=0;
    
    for(j=0;j<ntr;j++){
      
      if((i == 1 && j == idx[i-1]) 
	 ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
      
      if(charges[j] > qmax) {
	qmax = charges[j];
	jmax=j;
      }       
    } 
    
    if(qmax > 0) {
      idx[i]=jmax;
      jch[i]=charges[jmax]; 
      jtr[i]=tracks[jmax]; 
    }
    
  } 
  
  for(i=0;i<3;i++){
    if (jtr[i] == -2) {
         charges[i]=0;
         tracks[i]=0;
    } else {
         charges[i]=jch[i];
         tracks[i]=jtr[i];
    }
  }
}
Commit	Line	Data
66f93042	1	/**************************************************************************
	2	* Copyright(c) 1998-2000, 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	/*
	17	$Log$
116cbefd	18	Revision 1.9 2002/03/13 07:04:11 jchudoba
	19	Connect only MUON branches when reading the event to speed up digitisation.
	20
4b1670dc	21	Revision 1.8 2002/02/22 12:14:21 morsch
	22	Validate pad hit before digitization.
	23
1e2fa2a2	24	Revision 1.7 2001/11/22 11:15:41 jchudoba
	25	Proper deletion of arrays (thanks to Rene Brun)
	26
3ec837bb	27	Revision 1.6 2001/11/02 12:55:45 jchudoba
	28	cleanup of the code, add const to Get methods
	29
c7955c78	30	Revision 1.5 2001/10/31 16:35:07 jchudoba
	31	some functionality move to AliMUONTransientDigit class
	32
8f36c696	33	Revision 1.4 2001/03/20 13:36:11 egangler
	34	TFile memory leak and "too many files open" problem solved (?):
	35	the backgroud file is open only once forever, and not once for each event.
	36
77ccf5b0	37	Revision 1.3 2001/03/05 23:57:44 morsch
	38	Writing of digit tree moved to macro.
	39
7868442b	40	Revision 1.2 2001/03/05 08:40:25 morsch
	41	Method SortTracks(..) imported from AliMUON.
	42
a3f472a9	43	Revision 1.1 2001/02/02 14:11:53 morsch
	44	AliMUONMerger prototype to be called by the merge manager.
	45
66f93042	46	*/
66f93042	47
66f93042	48	#include <TDirectory.h>
116cbefd	49	#include <TFile.h>
	50	#include <TObjArray.h>
	51	#include <TPDGCode.h>
	52	#include <TTree.h>
66f93042	53
66f93042	54	#include "AliHitMap.h"
66f93042	55	#include "AliMUON.h"
116cbefd	56	#include "AliMUONChamber.h"
	57	#include "AliMUONConstants.h"
	58	#include "AliMUONDigit.h"
66f93042	59	#include "AliMUONHit.h"
116cbefd	60	#include "AliMUONHitMapA1.h"
116cbefd	61	#include "AliMUONMerger.h"
66f93042	62	#include "AliMUONPadHit.h"
66f93042	63	#include "AliMUONTransientDigit.h"
66f93042	64	#include "AliRun.h"
66f93042	65
	66	ClassImp(AliMUONMerger)
	67
	68	//___________________________________________
	69	AliMUONMerger::AliMUONMerger()
	70	{
	71	// Default constructor
	72	fEvNrSig = 0;
	73	fEvNrBgr = 0;
	74	fMerge = kDigitize;
	75	fFnBgr = 0;
	76	fHitMap = 0;
	77	fList = 0;
	78	fTrH1 = 0;
	79	fHitsBgr = 0;
	80	fPadHitsBgr = 0;
	81	fHitMap = 0;
	82	fList = 0;
77ccf5b0	83	fBgrFile = 0;
66f93042	84	}
	85
	86	//------------------------------------------------------------------------
	87	AliMUONMerger::~AliMUONMerger()
	88	{
	89	// Destructor
	90	if (fTrH1) delete fTrH1;
	91	if (fHitsBgr) delete fHitsBgr;
	92	if (fPadHitsBgr) delete fPadHitsBgr;
3ec837bb	93	if (fHitMap) delete [] fHitMap;
66f93042	94	if (fList) delete fList;
77ccf5b0	95	if (fBgrFile) delete fBgrFile;
66f93042	96	}
	97
	98	//------------------------------------------------------------------------
c7955c78	99	Bool_t AliMUONMerger::Exists(const AliMUONPadHit *padhit) const
66f93042	100	{
3ec837bb	101	// test if the given padhit was already fired
66f93042	102	return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
	103	}
	104
	105	//------------------------------------------------------------------------
8f36c696	106	void AliMUONMerger::Update(AliMUONPadHit *padhit)
66f93042	107	{
3ec837bb	108	// add new contribution to the fired padhit
8f36c696	109	AliMUONTransientDigit *pdigit =
	110	static_cast<AliMUONTransientDigit*>(
	111	fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
66f93042	112
66f93042	113	// update charge
66f93042	114	//
8f36c696	115	Int_t iqpad = padhit->QPad(); // charge per pad
	116	pdigit->AddSignal(iqpad);
	117	pdigit->AddPhysicsSignal(iqpad);
	118
66f93042	119	// update list of tracks
66f93042	120	//
8f36c696	121	Int_t track, charge;
	122	if (fSignal) {
	123	track = fTrack;
	124	charge = iqpad;
66f93042	125	} else {
8f36c696	126	track = kBgTag;
8f36c696	127	charge = kBgTag;
66f93042	128	}
8f36c696	129	pdigit->UpdateTrackList(track,charge);
66f93042	130	}
	131
	132	//------------------------------------------------------------------------
8f36c696	133	void AliMUONMerger::CreateNew(AliMUONPadHit *padhit)
66f93042	134	{
3ec837bb	135	// add new transient digit to the list, update hit map
66f93042	136	fList->AddAtAndExpand(
66f93042	137	new AliMUONTransientDigit(fNch,fDigits),fCounter);
8f36c696	138	fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
8f36c696	139	AliMUONTransientDigit* pdigit =
c7955c78	140	static_cast<AliMUONTransientDigit*>(fList->Last());
66f93042	141	// list of tracks
8f36c696	142	Int_t track, charge;
	143	if (fSignal) {
	144	track = fTrack;
	145	charge = padhit->QPad();
	146	} else {
	147	track = kBgTag;
	148	charge = kBgTag;
	149	}
	150	pdigit->AddToTrackList(track,charge);
	151	fCounter++;
66f93042	152	}
	153
	154
	155	//------------------------------------------------------------------------
	156	void AliMUONMerger::Init()
	157	{
	158	// Initialisation
77ccf5b0	159	// open only once the background file !!
77ccf5b0	160	if (fMerge && !fBgrFile) fBgrFile = InitBgr();
66f93042	161	}
	162
	163
	164
	165	//------------------------------------------------------------------------
	166	TFile* AliMUONMerger::InitBgr()
	167	{
	168	// Initialise background event
	169	TFile *file = new TFile(fFnBgr);
	170	// add error checking later
	171	printf("\n AliMUONMerger has opened %s file with background event \n", fFnBgr);
	172	fHitsBgr = new TClonesArray("AliMUONHit",1000);
	173	fPadHitsBgr = new TClonesArray("AliMUONPadHit",1000);
	174	return file;
	175	}
	176
	177	//------------------------------------------------------------------------
	178	void AliMUONMerger::Digitise()
	179	{
	180
	181	// keep galice.root for signal and name differently the file for
	182	// background when add! otherwise the track info for signal will be lost !
	183
	184	AliMUONChamber* iChamber;
	185	AliSegmentation* segmentation;
	186
	187	fList = new TObjArray;
66f93042	188
	189	AliMUON pMUON = (AliMUON ) gAlice->GetModule("MUON");
	190	fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
66f93042	191	if (fMerge ) {
66f93042	192	fBgrFile->cd();
66f93042	193	//
66f93042	194	// Get Hits Tree header from file
77ccf5b0	195	//if(fHitsBgr) fHitsBgr->Clear(); // Useless because line 327
77ccf5b0	196	//if(fPadHitsBgr) fPadHitsBgr->Clear(); // Useless because line 328
66f93042	197	if(fTrH1) delete fTrH1;
	198	fTrH1 = 0;
	199
	200	char treeName[20];
	201	sprintf(treeName,"TreeH%d",fEvNrBgr);
	202	fTrH1 = (TTree*)gDirectory->Get(treeName);
	203	if (!fTrH1) {
	204	printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
	205	}
	206	//
	207	// Set branch addresses
	208	TBranch *branch;
	209	char branchname[20];
	210	sprintf(branchname,"%s",pMUON->GetName());
	211	if (fTrH1 && fHitsBgr) {
	212	branch = fTrH1->GetBranch(branchname);
	213	if (branch) branch->SetAddress(&fHitsBgr);
	214	}
	215	if (fTrH1 && fPadHitsBgr) {
	216	branch = fTrH1->GetBranch("MUONCluster");
	217	if (branch) branch->SetAddress(&fPadHitsBgr);
	218	}
	219	}
	220	//
	221	// loop over cathodes
	222	//
8f36c696	223	fSignal = kTRUE;
66f93042	224	for (int icat = 0; icat < 2; icat++) {
66f93042	225	fCounter = 0;
66f93042	226	for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
	227	iChamber = &(pMUON->Chamber(i));
	228	if (iChamber->Nsec() == 1 && icat == 1) {
	229	continue;
	230	} else {
	231	segmentation = iChamber->SegmentationModel(icat+1);
	232	}
	233	fHitMap[i] = new AliMUONHitMapA1(segmentation, fList);
66f93042	234	}
	235
	236	//
	237	// Loop over tracks
	238	//
	239
	240	TTree *treeH = gAlice->TreeH();
	241	Int_t ntracks = (Int_t) treeH->GetEntries();
4b1670dc	242	treeH->SetBranchStatus("*",0); // switch off all branches
4b1670dc	243	treeH->SetBranchStatus("MUON*",1); // switch on only MUON
66f93042	244
	245	for (fTrack = 0; fTrack < ntracks; fTrack++) {
	246	gAlice->ResetHits();
4b1670dc	247	treeH->GetEntry(fTrack,0);
66f93042	248	//
	249	// Loop over hits
	250	for(AliMUONHit* mHit = (AliMUONHit*)pMUON->FirstHit(-1);
	251	mHit;
	252	mHit = (AliMUONHit*)pMUON->NextHit())
	253	{
	254	fNch = mHit->Chamber()-1; // chamber number
	255	if (fNch > AliMUONConstants::NCh()-1) continue;
	256	iChamber = &(pMUON->Chamber(fNch));
66f93042	257
	258	//
	259	// Loop over pad hits
	260	for (AliMUONPadHit* mPad =
	261	(AliMUONPadHit*)pMUON->FirstPad(mHit,pMUON->PadHits());
	262	mPad;
	263	mPad = (AliMUONPadHit*)pMUON->NextPad(pMUON->PadHits()))
	264	{
	265	Int_t cathode = mPad->Cathode(); // cathode number
66f93042	266	if (cathode != (icat+1)) continue;
c7955c78	267	Int_t iqpad = Int_t(mPad->QPad()); // charge per pad
	268	// segmentation = iChamber->SegmentationModel(cathode);
	269	fDigits[0] = mPad->PadX();
	270	fDigits[1] = mPad->PadY();
1e2fa2a2	271	if (!(fHitMap[fNch]->ValidateHit(fDigits[0], fDigits[1]))) continue;
66f93042	272	fDigits[2] = icat;
	273	fDigits[3] = iqpad;
	274	fDigits[4] = iqpad;
	275	if (mHit->Particle() == kMuonPlus \|\|
	276	mHit->Particle() == kMuonMinus) {
	277	fDigits[5] = mPad->HitNumber();
	278	} else fDigits[5] = -1;
	279
	280	// build the list of fired pads and update the info
	281
	282	if (!Exists(mPad)) {
	283	CreateNew(mPad);
	284	} else {
	285	Update(mPad);
	286	} // end if pdigit
	287	} //end loop over clusters
	288	} // hit loop
	289	} // track loop
	290
	291	// open the file with background
	292
	293	if (fMerge) {
8f36c696	294	fSignal = kFALSE;
66f93042	295	ntracks = (Int_t)fTrH1->GetEntries();
	296	//
	297	// Loop over tracks
	298	//
	299	for (fTrack = 0; fTrack < ntracks; fTrack++) {
	300
	301	if (fHitsBgr) fHitsBgr->Clear();
	302	if (fPadHitsBgr) fPadHitsBgr->Clear();
	303
	304	fTrH1->GetEvent(fTrack);
	305	//
	306	// Loop over hits
	307	AliMUONHit* mHit;
c7955c78	308	for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
66f93042	309	{
	310	mHit = (AliMUONHit) (fHitsBgr)[i];
	311	fNch = mHit->Chamber()-1; // chamber number
	312	iChamber = &(pMUON->Chamber(fNch));
66f93042	313	//
	314	// Loop over pad hits
	315	for (AliMUONPadHit* mPad =
	316	(AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHitsBgr);
	317	mPad;
	318	mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHitsBgr))
	319	{
	320	Int_t cathode = mPad->Cathode(); // cathode number
	321	Int_t ipx = mPad->PadX(); // pad number on X
	322	Int_t ipy = mPad->PadY(); // pad number on Y
	323	Int_t iqpad = Int_t(mPad->QPad()); // charge per pad
1e2fa2a2	324	if (!(fHitMap[fNch]->ValidateHit(ipx, ipy))) continue;
66f93042	325
66f93042	326	if (cathode != (icat+1)) continue;
66f93042	327	fDigits[0] = ipx;
	328	fDigits[1] = ipy;
	329	fDigits[2] = icat;
	330	fDigits[3] = iqpad;
	331	fDigits[4] = 0;
	332	fDigits[5] = -1;
	333
	334	// build the list of fired pads and update the info
	335	if (!Exists(mPad)) {
	336	CreateNew(mPad);
	337	} else {
	338	Update(mPad);
	339	} // end if !Exists
	340	} //end loop over clusters
	341	} // hit loop
	342	} // track loop
	343
c7955c78	344	TTree *treeK = gAlice->TreeK();
66f93042	345	TFile *file = NULL;
66f93042	346
c7955c78	347	if (treeK) file = treeK->GetCurrentFile();
66f93042	348	file->cd();
66f93042	349	} // if fMerge
66f93042	350
8f36c696	351	Int_t tracks[kMAXTRACKS];
8f36c696	352	Int_t charges[kMAXTRACKS];
66f93042	353	Int_t nentries = fList->GetEntriesFast();
	354
	355	for (Int_t nent = 0; nent < nentries; nent++) {
	356	AliMUONTransientDigit address = (AliMUONTransientDigit)fList->At(nent);
	357	if (address == 0) continue;
	358	Int_t ich = address->Chamber();
	359	Int_t q = address->Signal();
	360	iChamber = &(pMUON->Chamber(ich));
	361	//
	362	// Digit Response (noise, threshold, saturation, ...)
	363	AliMUONResponse * response = iChamber->ResponseModel();
	364	q = response->DigitResponse(q);
	365
	366	if (!q) continue;
	367
	368	fDigits[0] = address->PadX();
	369	fDigits[1] = address->PadY();
	370	fDigits[2] = address->Cathode();
	371	fDigits[3] = q;
	372	fDigits[4] = address->Physics();
	373	fDigits[5] = address->Hit();
	374
8f36c696	375	Int_t nptracks = address->GetNTracks();
66f93042	376
8f36c696	377	if (nptracks > kMAXTRACKS) {
	378	printf("\n Attention - nptracks > kMAXTRACKS %d \n", nptracks);
	379	nptracks = kMAXTRACKS;
66f93042	380	}
	381	if (nptracks > 2) {
	382	printf("Attention - nptracks > 2 %d \n",nptracks);
	383	printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
	384	}
	385	for (Int_t tr = 0; tr < nptracks; tr++) {
8f36c696	386	tracks[tr] = address->GetTrack(tr);
8f36c696	387	charges[tr] = address->GetCharge(tr);
66f93042	388	} //end loop over list of tracks for one pad
	389	// Sort list of tracks according to charge
	390	if (nptracks > 1) {
a3f472a9	391	SortTracks(tracks,charges,nptracks);
66f93042	392	}
8f36c696	393	if (nptracks < kMAXTRACKS ) {
8f36c696	394	for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
66f93042	395	tracks[i] = 0;
	396	charges[i] = 0;
	397	}
	398	}
	399
	400	// fill digits
	401	pMUON->AddDigits(ich,tracks,charges,fDigits);
66f93042	402	}
	403	gAlice->TreeD()->Fill();
	404	pMUON->ResetDigits();
	405	fList->Delete();
	406
	407
	408	for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
	409	if (fHitMap[ii]) {
c7955c78	410	delete fHitMap[ii];
66f93042	411	fHitMap[ii] = 0;
	412	}
	413	}
66f93042	414	} //end loop over cathodes
66f93042	415	delete [] fHitMap;
66f93042	416	delete fList;
66f93042	417
77ccf5b0	418	// no need to delete ... and it makes a crash also
	419	// if (fHitsBgr) fHitsBgr->Delete();
	420	// if (fPadHitsBgr) fPadHitsBgr->Delete();
66f93042	421	// gObjectTable->Print();
	422	}
	423
	424
	425
a3f472a9	426	void AliMUONMerger::SortTracks(Int_t tracks,Int_t charges,Int_t ntr)
	427	{
	428	//
	429	// Sort the list of tracks contributing to a given digit
	430	// Only the 3 most significant tracks are acctually sorted
	431	//
	432
	433	//
	434	// Loop over signals, only 3 times
	435	//
	436
	437	Int_t qmax;
	438	Int_t jmax;
	439	Int_t idx[3] = {-2,-2,-2};
	440	Int_t jch[3] = {-2,-2,-2};
	441	Int_t jtr[3] = {-2,-2,-2};
	442	Int_t i,j,imax;
	443
	444	if (ntr<3) imax=ntr;
	445	else imax=3;
	446	for(i=0;i<imax;i++){
	447	qmax=0;
	448	jmax=0;
	449
	450	for(j=0;j<ntr;j++){
	451
	452	if((i == 1 && j == idx[i-1])
	453	\|\|(i == 2 && (j == idx[i-1] \|\| j == idx[i-2]))) continue;
	454
	455	if(charges[j] > qmax) {
	456	qmax = charges[j];
	457	jmax=j;
	458	}
	459	}
	460
	461	if(qmax > 0) {
	462	idx[i]=jmax;
	463	jch[i]=charges[jmax];
	464	jtr[i]=tracks[jmax];
	465	}
	466
	467	}
	468
	469	for(i=0;i<3;i++){
	470	if (jtr[i] == -2) {
	471	charges[i]=0;
	472	tracks[i]=0;
	473	} else {
	474	charges[i]=jch[i];
	475	tracks[i]=jtr[i];
	476	}
	477	}
	478	}
	479
66f93042	480
66f93042	481