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

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


// #include "AliMerger.h"
// #include "AliMergable.h"
#include "AliMUONMerger.h"
#include "AliMUONConstants.h"
#include "AliMUONChamber.h"
#include "AliHitMap.h"
#include "AliMUONHitMapA1.h"
#include "AliMUON.h"
#include "AliMUONHit.h"
#include "AliMUONPadHit.h"
#include "AliMUONDigit.h"
#include "AliMUONTransientDigit.h"
#include "AliRun.h"
#include "AliPDG.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;
    fTrList     = 0;
    fAddress    = 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 (fTrList)     delete fTrList;
    if (fAddress)    delete fAddress; 
}

//------------------------------------------------------------------------
Bool_t AliMUONMerger::Exists(const AliMUONPadHit *mergable)
{
    AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;
    return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
}

//------------------------------------------------------------------------
void AliMUONMerger::Update(AliMUONPadHit *mergable)
{
    AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;    
    AliMUONTransientDigit* pdigit;
    Int_t ipx      = padhit->PadX();        // pad number on X
    Int_t ipy      = padhit->PadY();        // pad number on Y
    Int_t iqpad    = Int_t(padhit->QPad()); // charge per pad

    pdigit = (AliMUONTransientDigit*) fHitMap[fNch]->GetHit(ipx, ipy);
    // update charge
    //
    (*pdigit).AddSignal(iqpad);
    (*pdigit).AddPhysicsSignal(iqpad);		
    // update list of tracks
    //
    TObjArray* fTrList = (TObjArray*)pdigit->TrackList();
    Int_t lastEntry = fTrList->GetLast();
    TVector *pTrack = (TVector*)fTrList->At(lastEntry);
    TVector &ptrk   = *pTrack;
    TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
    Int_t lastTrack = Int_t(ptrk(0));

    if (trinfo(0) != kBgTag) {
	if (lastTrack == fTrack) {
	    Int_t lastCharge = Int_t(ptrk(1));
	    lastCharge += iqpad;
	    fTrList->RemoveAt(lastEntry);
	    trinfo(1) = lastCharge;
	    fTrList->AddAt(&trinfo,lastEntry);
	} else {
	    fTrList->Add(&trinfo);
	}
    } else {
	if (lastTrack != -1) fTrList->Add(&trinfo);
    }
}

//------------------------------------------------------------------------
void AliMUONMerger::CreateNew(AliMUONPadHit *mergable)
{
    AliMUONPadHit *padhit = (AliMUONPadHit*) mergable;    
    AliMUONTransientDigit* pdigit;

    Int_t ipx      = padhit->PadX();       // pad number on X
    Int_t ipy      = padhit->PadY();       // pad number on Y
    fList->AddAtAndExpand(
	new AliMUONTransientDigit(fNch,fDigits),fCounter);
    fHitMap[fNch]->SetHit(ipx, ipy, fCounter);
    fCounter++;
    pdigit = (AliMUONTransientDigit*)fList->At(fList->GetLast());
    // list of tracks
    TObjArray *fTrList = (TObjArray*)pdigit->TrackList();
    TVector &trinfo    = *((TVector*) (*fAddress)[fCountadr-1]);
    fTrList->Add(&trinfo);
}


//------------------------------------------------------------------------
void AliMUONMerger::Init()
{
// Initialisation
    
    if (fMerge) 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;
    if(!fAddress) fAddress = new TClonesArray("TVector",10000);

    AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
    fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
    for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {fHitMap[i] = 0;}
    if (fMerge ) {
	fBgrFile->cd();
	fBgrFile->ls();
        //
	// Get Hits Tree header from file
	if(fHitsBgr) fHitsBgr->Clear();
	if(fPadHitsBgr) fPadHitsBgr->Clear();
	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
    //
    AliHitMap* hm;
    fCountadr = 0;
    for (int icat = 0; icat < 2; icat++) { 
	fCounter = 0;
	Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
	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);
	    nmuon[i] = 0;
	}

//
//   Loop over tracks
//

	TTree *treeH  = gAlice->TreeH();
	Int_t ntracks = (Int_t) treeH->GetEntries();
	Int_t jj;

	Float_t ** xhit = new Float_t * [AliMUONConstants::NCh()];
	for (jj = 0; jj < AliMUONConstants::NCh(); jj++) xhit[jj] = new Float_t[2];
	Float_t ** yhit = new Float_t * [AliMUONConstants::NCh()];
	for (jj = 0; jj < AliMUONConstants::NCh(); jj++) yhit[jj] = new Float_t[2];

	for (fTrack = 0; fTrack < ntracks; fTrack++) {
	    gAlice->ResetHits();
	    treeH->GetEvent(fTrack);
//
//   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));
		/*
		if (fMerge) {
		    if (mHit->Particle() == kMuonPlus || 
			mHit->Particle() == kMuonMinus) {
			xhit[fNch][nmuon[fNch]] = mHit->X();
			yhit[fNch][nmuon[fNch]] = mHit->Y();
			nmuon[fNch]++;
			if (nmuon[fNch] > 2) printf("MUON: nmuon %d\n",nmuon[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
		    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 (cathode != (icat+1)) continue;

		    segmentation = iChamber->SegmentationModel(cathode);

		    new((*fAddress)[fCountadr++]) TVector(2);

		    TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
		    trinfo(0) = (Float_t)fTrack;
		    trinfo(1) = (Float_t)iqpad;

		    fDigits[0] = ipx;
		    fDigits[1] = ipy;
		    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) {
	    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 i = 0; i < fHitsBgr->GetEntriesFast(); ++i) 
		{	
		    mHit   = (AliMUONHit*) (*fHitsBgr)[i];
		    fNch   = mHit->Chamber()-1;  // chamber number
		    iChamber = &(pMUON->Chamber(fNch));
                    Float_t xbgr = mHit->X();
		    Float_t ybgr = mHit->Y();
		    Bool_t cond  = kFALSE;
		    /*
		    for (Int_t imuon = 0; imuon < nmuon[fNch]; imuon++) {
			Float_t dist = (xbgr-xhit[fNch][imuon])*(xbgr-xhit[fNch][imuon])
			    +(ybgr-yhit[fNch][imuon])*(ybgr-yhit[fNch][imuon]);
			if (dist < 100.) cond = kTRUE;
		    }
		    */
		    cond  = kTRUE;
//
// 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 (cathode != (icat+1)) continue;
			new((*fAddress)[fCountadr++]) TVector(2);
			TVector &trinfo = *((TVector*) (*fAddress)[fCountadr-1]);
			trinfo(0) = kBgTag;  // tag background
			trinfo(1) = kBgTag;
			
			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 *fAli = gAlice->TreeK();
            TFile *file = NULL;
	    
	    if (fAli) file = fAli->GetCurrentFile();
	    file->cd();
	} // if fMerge
	delete [] xhit;
	delete [] yhit;

	Int_t tracks[10];
	Int_t charges[10];
	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();
	    
	    TObjArray* fTrList = (TObjArray*)address->TrackList();
	    Int_t nptracks = fTrList->GetEntriesFast();

	    // this was changed to accomodate the real number of tracks

	    if (nptracks > 10) {
		printf("\n Attention - nptracks > 10 %d \n", nptracks);
		nptracks = 10;
	    }
	    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++) {
		TVector *ppP = (TVector*)fTrList->At(tr);
		if(!ppP ) printf("ppP - %p\n",ppP);
		TVector &pp  = *ppP;
		tracks[tr]   = Int_t(pp(0));
		charges[tr]  = Int_t(pp(1));
	    }      //end loop over list of tracks for one pad
            // Sort list of tracks according to charge
	    if (nptracks > 1) {
		pMUON->SortTracks(tracks,charges,nptracks);
	    }
	    if (nptracks < 10 ) {
		for (Int_t i = nptracks; i < 10; i++) {
		    tracks[i]  = 0;
		    charges[i] = 0;
		}
	    }
	    
	    // fill digits
	    pMUON->AddDigits(ich,tracks,charges,fDigits);
	    // delete fTrList;
	}
	gAlice->TreeD()->Fill();
	pMUON->ResetDigits();
	fList->Delete();

	
	for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
	    if (fHitMap[ii]) {
		hm=fHitMap[ii];
		delete hm;
		fHitMap[ii] = 0;
	    }
	}
	delete [] nmuon;    
    } //end loop over cathodes
    delete [] fHitMap;
    char hname[30];
    sprintf(hname,"TreeD%d",fEvNrSig);
    gAlice->TreeD()->Write(hname);
    // reset tree
    gAlice->TreeD()->Reset();
    delete fList;
    
    if (fAddress)    fAddress->Delete();
    if (fHitsBgr)    fHitsBgr->Delete();
    if (fPadHitsBgr) fPadHitsBgr->Delete();
    // gObjectTable->Print();
}
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$
	18	*/
	19
	20	#include <TTree.h>
	21	#include <TVector.h>
	22	#include <TObjArray.h>
	23	#include <TFile.h>
	24	#include <TDirectory.h>
	25
	26
	27	// #include "AliMerger.h"
	28	// #include "AliMergable.h"
	29	#include "AliMUONMerger.h"
	30	#include "AliMUONConstants.h"
	31	#include "AliMUONChamber.h"
	32	#include "AliHitMap.h"
	33	#include "AliMUONHitMapA1.h"
	34	#include "AliMUON.h"
	35	#include "AliMUONHit.h"
	36	#include "AliMUONPadHit.h"
	37	#include "AliMUONDigit.h"
	38	#include "AliMUONTransientDigit.h"
	39	#include "AliRun.h"
	40	#include "AliPDG.h"
	41
	42	ClassImp(AliMUONMerger)
	43
	44	//___________________________________________
	45	AliMUONMerger::AliMUONMerger()
	46	{
	47	// Default constructor
	48	fEvNrSig = 0;
	49	fEvNrBgr = 0;
	50	fMerge = kDigitize;
	51	fFnBgr = 0;
	52	fHitMap = 0;
	53	fList = 0;
	54	fTrH1 = 0;
	55	fHitsBgr = 0;
	56	fPadHitsBgr = 0;
	57	fHitMap = 0;
	58	fList = 0;
	59	fTrList = 0;
	60	fAddress = 0;
	61	}
	62
	63	//------------------------------------------------------------------------
	64	AliMUONMerger::~AliMUONMerger()
65	{
66	// Destructor
67	if (fTrH1) delete fTrH1;
68	if (fHitsBgr) delete fHitsBgr;
69	if (fPadHitsBgr) delete fPadHitsBgr;
70	if (fHitMap) delete fHitMap;
71	if (fList) delete fList;
72	if (fTrList) delete fTrList;
73	if (fAddress) delete fAddress;
74	}
75
76	//------------------------------------------------------------------------
77	Bool_t AliMUONMerger::Exists(const AliMUONPadHit *mergable)
78	{
79	AliMUONPadHit padhit = (AliMUONPadHit) mergable;
80	return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
81	}
82
83	//------------------------------------------------------------------------
84	void AliMUONMerger::Update(AliMUONPadHit *mergable)
85	{
86	AliMUONPadHit padhit = (AliMUONPadHit) mergable;
87	AliMUONTransientDigit* pdigit;
88	Int_t ipx = padhit->PadX(); // pad number on X
89	Int_t ipy = padhit->PadY(); // pad number on Y
90	Int_t iqpad = Int_t(padhit->QPad()); // charge per pad
91
92	pdigit = (AliMUONTransientDigit*) fHitMap[fNch]->GetHit(ipx, ipy);
93	// update charge
94	//
95	(*pdigit).AddSignal(iqpad);
96	(*pdigit).AddPhysicsSignal(iqpad);
97	// update list of tracks
98	//
99	TObjArray* fTrList = (TObjArray*)pdigit->TrackList();
100	Int_t lastEntry = fTrList->GetLast();
101	TVector pTrack = (TVector)fTrList->At(lastEntry);
102	TVector &ptrk = *pTrack;
103	TVector &trinfo = ((TVector) (*fAddress)[fCountadr-1]);
104	Int_t lastTrack = Int_t(ptrk(0));
105
106	if (trinfo(0) != kBgTag) {
107	if (lastTrack == fTrack) {
108	Int_t lastCharge = Int_t(ptrk(1));
109	lastCharge += iqpad;
110	fTrList->RemoveAt(lastEntry);
111	trinfo(1) = lastCharge;
112	fTrList->AddAt(&trinfo,lastEntry);
113	} else {
114	fTrList->Add(&trinfo);
115	}
116	} else {
117	if (lastTrack != -1) fTrList->Add(&trinfo);
118	}
119	}
120
121	//------------------------------------------------------------------------
122	void AliMUONMerger::CreateNew(AliMUONPadHit *mergable)
123	{
124	AliMUONPadHit padhit = (AliMUONPadHit) mergable;
125	AliMUONTransientDigit* pdigit;
126
127	Int_t ipx = padhit->PadX(); // pad number on X
128	Int_t ipy = padhit->PadY(); // pad number on Y
129	fList->AddAtAndExpand(
130	new AliMUONTransientDigit(fNch,fDigits),fCounter);
131	fHitMap[fNch]->SetHit(ipx, ipy, fCounter);
132	fCounter++;
133	pdigit = (AliMUONTransientDigit*)fList->At(fList->GetLast());
134	// list of tracks
135	TObjArray fTrList = (TObjArray)pdigit->TrackList();
136	TVector &trinfo = ((TVector) (*fAddress)[fCountadr-1]);
137	fTrList->Add(&trinfo);
138	}
139
140
141	//------------------------------------------------------------------------
142	void AliMUONMerger::Init()
143	{
144	// Initialisation
145
146	if (fMerge) fBgrFile = InitBgr();
147	}
148
149
150
151	//------------------------------------------------------------------------
152	TFile* AliMUONMerger::InitBgr()
153	{
154	// Initialise background event
155	TFile *file = new TFile(fFnBgr);
156	// add error checking later
157	printf("\n AliMUONMerger has opened %s file with background event \n", fFnBgr);
158	fHitsBgr = new TClonesArray("AliMUONHit",1000);
159	fPadHitsBgr = new TClonesArray("AliMUONPadHit",1000);
160	return file;
161	}
162
163	//------------------------------------------------------------------------
164	void AliMUONMerger::Digitise()
165	{
166
167	// keep galice.root for signal and name differently the file for
168	// background when add! otherwise the track info for signal will be lost !
169
170	AliMUONChamber* iChamber;
171	AliSegmentation* segmentation;
172
173	fList = new TObjArray;
174	if(!fAddress) fAddress = new TClonesArray("TVector",10000);
175
176	AliMUON pMUON = (AliMUON ) gAlice->GetModule("MUON");
177	fHitMap= new AliHitMap* [AliMUONConstants::NCh()];
178	for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {fHitMap[i] = 0;}
179	if (fMerge ) {
180	fBgrFile->cd();
181	fBgrFile->ls();
182	//
183	// Get Hits Tree header from file
184	if(fHitsBgr) fHitsBgr->Clear();
185	if(fPadHitsBgr) fPadHitsBgr->Clear();
186	if(fTrH1) delete fTrH1;
187	fTrH1 = 0;
188
189	char treeName[20];
190	sprintf(treeName,"TreeH%d",fEvNrBgr);
191	fTrH1 = (TTree*)gDirectory->Get(treeName);
192	if (!fTrH1) {
193	printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
194	}
195	//
196	// Set branch addresses
197	TBranch *branch;
198	char branchname[20];
199	sprintf(branchname,"%s",pMUON->GetName());
200	if (fTrH1 && fHitsBgr) {
201	branch = fTrH1->GetBranch(branchname);
202	if (branch) branch->SetAddress(&fHitsBgr);
203	}
204	if (fTrH1 && fPadHitsBgr) {
205	branch = fTrH1->GetBranch("MUONCluster");
206	if (branch) branch->SetAddress(&fPadHitsBgr);
207	}
208	}
209	//
210	// loop over cathodes
211	//
212	AliHitMap* hm;
213	fCountadr = 0;
214	for (int icat = 0; icat < 2; icat++) {
215	fCounter = 0;
216	Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
217	for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
218	iChamber = &(pMUON->Chamber(i));
219	if (iChamber->Nsec() == 1 && icat == 1) {
220	continue;
221	} else {
222	segmentation = iChamber->SegmentationModel(icat+1);
223	}
224	fHitMap[i] = new AliMUONHitMapA1(segmentation, fList);
225	nmuon[i] = 0;
226	}
227
228	//
229	// Loop over tracks
230	//
231
232	TTree *treeH = gAlice->TreeH();
233	Int_t ntracks = (Int_t) treeH->GetEntries();
234	Int_t jj;
235
236	Float_t ** xhit = new Float_t * [AliMUONConstants::NCh()];
237	for (jj = 0; jj < AliMUONConstants::NCh(); jj++) xhit[jj] = new Float_t[2];
238	Float_t ** yhit = new Float_t * [AliMUONConstants::NCh()];
239	for (jj = 0; jj < AliMUONConstants::NCh(); jj++) yhit[jj] = new Float_t[2];
240
241	for (fTrack = 0; fTrack < ntracks; fTrack++) {
242	gAlice->ResetHits();
243	treeH->GetEvent(fTrack);
244	//
245	// Loop over hits
246	for(AliMUONHit* mHit = (AliMUONHit*)pMUON->FirstHit(-1);
247	mHit;
248	mHit = (AliMUONHit*)pMUON->NextHit())
249	{
250	fNch = mHit->Chamber()-1; // chamber number
251	if (fNch > AliMUONConstants::NCh()-1) continue;
252	iChamber = &(pMUON->Chamber(fNch));
253	/*
254	if (fMerge) {
255	if (mHit->Particle() == kMuonPlus \|\|
256	mHit->Particle() == kMuonMinus) {
257	xhit[fNch][nmuon[fNch]] = mHit->X();
258	yhit[fNch][nmuon[fNch]] = mHit->Y();
259	nmuon[fNch]++;
260	if (nmuon[fNch] > 2) printf("MUON: nmuon %d\n",nmuon[fNch]);
261	}
262	}
263	*/
264
265	//
266	// Loop over pad hits
267	for (AliMUONPadHit* mPad =
268	(AliMUONPadHit*)pMUON->FirstPad(mHit,pMUON->PadHits());
269	mPad;
270	mPad = (AliMUONPadHit*)pMUON->NextPad(pMUON->PadHits()))
271	{
272	Int_t cathode = mPad->Cathode(); // cathode number
273	Int_t ipx = mPad->PadX(); // pad number on X
274	Int_t ipy = mPad->PadY(); // pad number on Y
275	Int_t iqpad = Int_t(mPad->QPad()); // charge per pad
276	if (cathode != (icat+1)) continue;
277
278	segmentation = iChamber->SegmentationModel(cathode);
279
280	new((*fAddress)[fCountadr++]) TVector(2);
281
282	TVector &trinfo = ((TVector) (*fAddress)[fCountadr-1]);
283	trinfo(0) = (Float_t)fTrack;
284	trinfo(1) = (Float_t)iqpad;
285
286	fDigits[0] = ipx;
287	fDigits[1] = ipy;
288	fDigits[2] = icat;
289	fDigits[3] = iqpad;
290	fDigits[4] = iqpad;
291	if (mHit->Particle() == kMuonPlus \|\|
292	mHit->Particle() == kMuonMinus) {
293	fDigits[5] = mPad->HitNumber();
294	} else fDigits[5] = -1;
295
296	// build the list of fired pads and update the info
297
298	if (!Exists(mPad)) {
299	CreateNew(mPad);
300	} else {
301	Update(mPad);
302	} // end if pdigit
303	} //end loop over clusters
304	} // hit loop
305	} // track loop
306
307	// open the file with background
308
309	if (fMerge) {
310	ntracks = (Int_t)fTrH1->GetEntries();
311	//
312	// Loop over tracks
313	//
314	for (fTrack = 0; fTrack < ntracks; fTrack++) {
315
316	if (fHitsBgr) fHitsBgr->Clear();
317	if (fPadHitsBgr) fPadHitsBgr->Clear();
318
319	fTrH1->GetEvent(fTrack);
320	//
321	// Loop over hits
322	AliMUONHit* mHit;
323	for(int i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
324	{
325	mHit = (AliMUONHit) (fHitsBgr)[i];
326	fNch = mHit->Chamber()-1; // chamber number
327	iChamber = &(pMUON->Chamber(fNch));
328	Float_t xbgr = mHit->X();
329	Float_t ybgr = mHit->Y();
330	Bool_t cond = kFALSE;
331	/*
332	for (Int_t imuon = 0; imuon < nmuon[fNch]; imuon++) {
333	Float_t dist = (xbgr-xhit[fNch][imuon])*(xbgr-xhit[fNch][imuon])
334	+(ybgr-yhit[fNch][imuon])*(ybgr-yhit[fNch][imuon]);
335	if (dist < 100.) cond = kTRUE;
336	}
337	*/
338	cond = kTRUE;
339	//
340	// Loop over pad hits
341	for (AliMUONPadHit* mPad =
342	(AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHitsBgr);
343	mPad;
344	mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHitsBgr))
345	{
346	Int_t cathode = mPad->Cathode(); // cathode number
347	Int_t ipx = mPad->PadX(); // pad number on X
348	Int_t ipy = mPad->PadY(); // pad number on Y
349	Int_t iqpad = Int_t(mPad->QPad()); // charge per pad
350
351	if (cathode != (icat+1)) continue;
352	new((*fAddress)[fCountadr++]) TVector(2);
353	TVector &trinfo = ((TVector) (*fAddress)[fCountadr-1]);
354	trinfo(0) = kBgTag; // tag background
355	trinfo(1) = kBgTag;
356
357	fDigits[0] = ipx;
358	fDigits[1] = ipy;
359	fDigits[2] = icat;
360	fDigits[3] = iqpad;
361	fDigits[4] = 0;
362	fDigits[5] = -1;
363
364	// build the list of fired pads and update the info
365	if (!Exists(mPad)) {
366	CreateNew(mPad);
367	} else {
368	Update(mPad);
369	} // end if !Exists
370	} //end loop over clusters
371	} // hit loop
372	} // track loop
373
374	TTree *fAli = gAlice->TreeK();
375	TFile *file = NULL;
376
377	if (fAli) file = fAli->GetCurrentFile();
378	file->cd();
379	} // if fMerge
380	delete [] xhit;
381	delete [] yhit;
382
383	Int_t tracks[10];
384	Int_t charges[10];
385	Int_t nentries = fList->GetEntriesFast();
386
387	for (Int_t nent = 0; nent < nentries; nent++) {
388	AliMUONTransientDigit address = (AliMUONTransientDigit)fList->At(nent);
389	if (address == 0) continue;
390	Int_t ich = address->Chamber();
391	Int_t q = address->Signal();
392	iChamber = &(pMUON->Chamber(ich));
393	//
394	// Digit Response (noise, threshold, saturation, ...)
395	AliMUONResponse * response = iChamber->ResponseModel();
396	q = response->DigitResponse(q);
397
398	if (!q) continue;
399
400	fDigits[0] = address->PadX();
401	fDigits[1] = address->PadY();
402	fDigits[2] = address->Cathode();
403	fDigits[3] = q;
404	fDigits[4] = address->Physics();
405	fDigits[5] = address->Hit();
406
407	TObjArray* fTrList = (TObjArray*)address->TrackList();
408	Int_t nptracks = fTrList->GetEntriesFast();
409
410	// this was changed to accomodate the real number of tracks
411
412	if (nptracks > 10) {
413	printf("\n Attention - nptracks > 10 %d \n", nptracks);
414	nptracks = 10;
415	}
416	if (nptracks > 2) {
417	printf("Attention - nptracks > 2 %d \n",nptracks);
418	printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
419	}
420	for (Int_t tr = 0; tr < nptracks; tr++) {
421	TVector ppP = (TVector)fTrList->At(tr);
422	if(!ppP ) printf("ppP - %p\n",ppP);
423	TVector &pp = *ppP;
424	tracks[tr] = Int_t(pp(0));
425	charges[tr] = Int_t(pp(1));
426	} //end loop over list of tracks for one pad
427	// Sort list of tracks according to charge
428	if (nptracks > 1) {
429	pMUON->SortTracks(tracks,charges,nptracks);
430	}
431	if (nptracks < 10 ) {
432	for (Int_t i = nptracks; i < 10; i++) {
433	tracks[i] = 0;
434	charges[i] = 0;
435	}
436	}
437
438	// fill digits
439	pMUON->AddDigits(ich,tracks,charges,fDigits);
440	// delete fTrList;
441	}
442	gAlice->TreeD()->Fill();
443	pMUON->ResetDigits();
444	fList->Delete();
445
446
447	for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
448	if (fHitMap[ii]) {
449	hm=fHitMap[ii];
450	delete hm;
451	fHitMap[ii] = 0;
452	}
453	}
454	delete [] nmuon;
455	} //end loop over cathodes
456	delete [] fHitMap;
457	char hname[30];
458	sprintf(hname,"TreeD%d",fEvNrSig);
459	gAlice->TreeD()->Write(hname);
460	// reset tree
461	gAlice->TreeD()->Reset();
462	delete fList;
463
464	if (fAddress) fAddress->Delete();
465	if (fHitsBgr) fHitsBgr->Delete();
466	if (fPadHitsBgr) fPadHitsBgr->Delete();
467	// gObjectTable->Print();
468	}
469
470
471
472
473