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



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

#include "AliMUON.h"
#include "AliMUONChamber.h"
#include "AliMUONConstants.h"
#include "AliMUONDigit.h"
#include "AliMUONSDigitizerv1.h"
#include "AliMUONHit.h"
#include "AliMUONHitMapA1.h"
#include "AliMUONPadHit.h"
#include "AliMUONTransientDigit.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliLoader.h"

ClassImp(AliMUONSDigitizerv1)

//___________________________________________
AliMUONSDigitizerv1::AliMUONSDigitizerv1() :
  fHitMap(0),
  fTDList(0),
  fTDCounter(0),
  fDebug(0),
  fMask(0),
  fSignal(0)
{
// Default ctor - don't use it
  if (GetDebug()>2) 
    cerr<<"AliMUONSDigitizerv1::AliMUONSDigitizerv1"
	<<"(AliRunDigitizer* manager) was processed"<<endl;
}


//------------------------------------------------------------------------
AliMUONSDigitizerv1::~AliMUONSDigitizerv1()
{
// Destructor
}

//------------------------------------------------------------------------
void AliMUONSDigitizerv1::AddTransientDigit(AliMUONTransientDigit * mTD)
{
  // Choosing the maping of the cathode plane of the chamber:
  Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
  fTDList->AddAtAndExpand(mTD, fTDCounter);
  fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
  fTDCounter++;
}

//------------------------------------------------------------------------
Bool_t AliMUONSDigitizerv1::ExistTransientDigit(AliMUONTransientDigit * mTD)
{
  // Choosing the maping of the cathode plane of the chamber:
  Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
  return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
}

//------------------------------------------------------------------------
Bool_t AliMUONSDigitizerv1::Init()
{

// Initialization 
  if (GetDebug()>2) Info("Init","AliMUONSDigitizerv1::Init() starts");

  //Loaders (We assume input0 to be the output too)
  AliRunLoader * runloader;    // Input   loader
  AliLoader    * gime; 
 
  // Getting runloader
  runloader    = AliRunLoader::GetRunLoader();
  if (runloader == 0x0) {
    Error("Init","RunLoader is not in input file 0");
    return kFALSE;  // RunDigitizer is not working.  
  }
  // Getting MUONloader
  gime        = runloader->GetLoader("MUONLoader");
  gime->LoadHits("READ");
  gime->LoadSDigits("RECREATE");
 
  return kTRUE;
}

//------------------------------------------------------------------------
void AliMUONSDigitizerv1::UpdateTransientDigit(Int_t track, AliMUONTransientDigit * mTD)
{
  // Choosing the maping of the cathode plane of the chamber:
  Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
  AliMUONTransientDigit *pdigit = 
    static_cast<AliMUONTransientDigit*>(fHitMap[iNchCpl]->GetHit(mTD->PadX(),mTD->PadY()));
  // update charge
  //
  Int_t iqpad  = mTD->Signal();        // charge per pad
  pdigit->AddSignal(iqpad);
  pdigit->AddPhysicsSignal(iqpad);		
  // update list of tracks
  //
  Int_t charge;   
  track=+ fMask;
  if (fSignal) charge = iqpad;
  //else         charge = kBgTag;
  else         charge = iqpad + fMask; 
	       
  pdigit->UpdateTrackList(track,charge);
}


//--------------------------------------------------------------------------
void AliMUONSDigitizerv1::MakeTransientDigit(Int_t track, Int_t iHit, AliMUONHit * mHit)
{
  AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
  if (!pMUON) { 
    cerr<<"AliMUONSDigitizerv1::MakeTransientDigit Error:"
	<<" module MUON not found in the input file"<<endl;
  }
  if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit starts"<<endl;
  Int_t   ichamber = mHit->Chamber()-1;
  AliMUONChamber & chamber = pMUON->Chamber(ichamber);
  Float_t xhit = mHit->X();
  Float_t yhit = mHit->Y();
  Float_t zhit = mHit->Z();
  Float_t eloss= mHit->Eloss(); 
  Float_t tof  = mHit->Age();
  // Variables for chamber response from AliMUONChamber::DisIntegration
  Float_t newdigit[6][500];  // Pad information
  Int_t nnew=0;              // Number of touched Pads per hit
  Int_t digits[6];
  
  //
  // Calls the charge disintegration method of the current chamber 
  if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit calling AliMUONChamber::DisIngtegration starts"<<endl;
  chamber.DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newdigit);
  // Creating a new TransientDigits from hit
  for(Int_t iTD=0; iTD<nnew; iTD++) {
    digits[0] = Int_t(newdigit[1][iTD]);  // Padx of the Digit
    digits[1] = Int_t(newdigit[2][iTD]);  // Pady of the Digit
    digits[2] = Int_t(newdigit[5][iTD]);  // Cathode plane
    digits[3] = Int_t(newdigit[3][iTD]);  // Induced charge in the Pad !!!! Int_t not correct
    if (fSignal) digits[4] = Int_t(newdigit[3][iTD]);
    else         digits[4] = 0;
    digits[5] = iHit+fMask;    // Hit number in the list
    if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit " <<
    cerr<<"AliMUONSDigitizerv1::MakeTransientDigit " <<
			"PadX "<< digits[0] << " " <<
			"PadY "<< digits[1] << " " <<
			"Plane " << digits[2] << " " <<
			"Charge " << digits[3] <<" " <<
			"newdigit " << newdigit[3][iTD] <<" " <<
			"Hit " << digits[5] << endl;
    // list of tracks
    Int_t charge;   
    track += fMask;
    if (fSignal) charge = digits[3];
    //else         charge = kBgTag;
    else         charge = digits[3] + fMask;
		 
    if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit Creating AliMUONTransientDigit"<<endl;
    AliMUONTransientDigit * mTD = new AliMUONTransientDigit(ichamber, digits);
    mTD->AddToTrackList(track,charge);
    if (!ExistTransientDigit(mTD)) {
      AddTransientDigit(mTD);
      if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit Adding TransientDigit"<<endl;
    }
    else {
      if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit updating TransientDigit"<<endl;
      UpdateTransientDigit(track, mTD);
      delete mTD;
    }
  }
}
//-----------------------------------------------------------------------
void AliMUONSDigitizerv1::Exec(Option_t* option)
{
  TString optionString = option;
  if (optionString.Data() == "deb") {
    Info("SDigitize","Called with option deb ");
    fDebug = 3;
  }

  AliMUONChamber*   chamber;
  AliSegmentation*  c1Segmentation; //Cathode plane c1 of the chamber
  AliSegmentation*  c2Segmentation; //Cathode place c2 of the chamber

  if (GetDebug()>2) Info("SDigitize","AliMUONSDigitizerv1::Exec starts");
  fTDList = new TObjArray;

  //Loaders (We assume input0 to be the output too)
  AliRunLoader * runloader;    // Input   loader
  AliLoader    * gime; 
 
  // Getting runloader
  runloader    = AliRunLoader::GetRunLoader();
  if (runloader == 0x0) {
    Error("SDigitize","RunLoader is not in input file 0");
    return;  // RunDigitizer is not working.  
  }
  // Getting MUONloader
  gime        = runloader->GetLoader("MUONLoader");
  if (gime->TreeH()==0x0) {
      if (GetDebug()>2) Info("SDigitize","TreeH is not loaded yet. Loading...");
     gime->LoadHits("READ");
       if (GetDebug()>2) Info("SDigitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
  }

  if (GetDebug()>2) Info("SDigitize","Loaders ready");

  if (runloader->GetAliRun() == 0x0) runloader->LoadgAlice();
  gAlice = runloader->GetAliRun();

  // Getting Module MUON  
  AliMUON *pMUON  = (AliMUON *) gAlice->GetDetector("MUON");
  if (!pMUON) {
    Error("SDigitize","Module MUON not found in the input file");
    return;
  }
  // Getting Muon data
  AliMUONData * muondata = pMUON->GetMUONData(); 
  muondata->SetLoader(gime);
  muondata->SetTreeAddress("H");

  Int_t currentevent = runloader->GetEventNumber();
  
  if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::Exec Event Number is "<<currentevent <<endl;
  if ( (currentevent<10)                                                 || 
       (Int_t(TMath::Log10(currentevent)) == TMath::Log10(currentevent) ) )
    cout <<"AliMUONSDigitizerv1::Exec Event Number is "<< currentevent <<endl;

  // New branch per chamber for MUON digit in the tree of digits
  if (gime->TreeS() == 0x0) {
    gime->MakeSDigitsContainer();
  }
  TTree* treeS = gime->TreeS();
  muondata->MakeBranch("S");
  muondata->SetTreeAddress("S");

  // Array of pointer of the AliMUONHitMapA1:
  //  two HitMaps per chamber, or one HitMap per cahtode plane
  fHitMap= new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];

  //Loop over chambers for the definition AliMUONHitMap
  for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
    chamber = &(pMUON->Chamber(i));
    c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
    fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
    c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
    fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
  }

// Loop to Sdigitize
    fSignal = kTRUE;
    
    // Setting the address 
    TTree *treeH = gime->TreeH();
    if (treeH == 0x0) {
      Error("SDigitize","Can not get TreeH from input ");
      Info("SDigitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
      return;
    }
    if (GetDebug()>2) {
      cerr<<"AliMUONSDigitizerv1::Exec treeH" << treeH <<endl;
    }
	
    if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::Exec Setting tree addresses"<<endl;

    //
    // Loop over tracks
    Int_t itrack;
    Int_t ntracks = (Int_t) treeH->GetEntries();
    for (itrack = 0; itrack < ntracks; itrack++) {
      if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::Exec itrack = "<<itrack<<endl;
      muondata->ResetHits();
      treeH->GetEvent(itrack);
      //
      //  Loop over hits
      Int_t ihit, ichamber;
      AliMUONHit* mHit;
      TClonesArray* hits = muondata->Hits();
      for(ihit = 0; ihit < hits->GetEntriesFast(); ihit++) {
	mHit = static_cast<AliMUONHit*>(hits->At(ihit));
	ichamber = mHit->Chamber()-1;  // chamber number
	if (ichamber > AliMUONConstants::NCh()-1) {
	  cerr<<"AliMUONSDigitizer: ERROR: "
	      <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
	      <<ichamber<<", "<< AliMUONConstants::NCh()<<endl;
	  return;
	}
	chamber = &(pMUON->Chamber(ichamber));
	//
	//Dumping Hit content:
	if (GetDebug()>2) {
	  cerr<<"AliMuonDigitizerv1::Exec ihit, ichamber, x, y, z, eloss " <<
	    ihit << " " << 
	    mHit->Chamber() << " " <<
	    mHit->X() << " " <<
	    mHit->Y() << " " <<
	    mHit->Z() << " " <<
	    mHit->Eloss() << " " << endl;
	}
	// 
	// Inititializing Correlation
	chamber->ChargeCorrelationInit();
	if (ichamber < AliMUONConstants::NTrackingCh()) {
	  // Tracking Chamber
	  // Initialize hit position (cursor) in the segmentation model 
	  chamber->SigGenInit(mHit->X(), mHit->Y(), mHit->Z());
	} else {
	  // Trigger Chamber 
	}
	MakeTransientDigit(itrack, ihit, mHit);
      } // hit loop
    } // track loop
    if (GetDebug()>2) cerr<<"AliMUONSDigitizer::Exec End of hits, track and file loops"<<endl;

    // Loop on cathodes
    Int_t icat;
    for(icat=0; icat<2; icat++) {
      //
      // Filling SDigit List
      Int_t tracks[kMAXTRACKS];
      Int_t charges[kMAXTRACKS];
      Int_t nentries = fTDList->GetEntriesFast();
      Int_t digits[6];
      for (Int_t nent = 0; nent < nentries; nent++) {
	AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
	if (address == 0) continue; 
	Int_t ich = address->Chamber();
	Int_t   q = address->Signal(); 
	
	if (!q) continue; // not mandatory ?
	
	digits[0] = address->PadX();
	digits[1] = address->PadY();
	digits[2] = address->Cathode()-1;
	digits[3] = q;
	digits[4] = address->Physics();
	digits[5] = address->Hit();
	
	Int_t nptracks = address->GetNTracks();
	
	if (nptracks > kMAXTRACKS) {
	  if (GetDebug() >0) {
	    cerr<<"AliMUONSDigitizer:Exec  nptracks > 10 "<<nptracks;
	    cerr<<"reset to max value "<<kMAXTRACKS<<endl;
	  }
	  nptracks = kMAXTRACKS;
	}
	if (nptracks > 2 && GetDebug() >2) {
	  cerr<<"AliMUONSDigitizer::Exec  nptracks > 2 "<<nptracks<<endl;
	  //	printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[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]  = -1;
	    charges[i] = 0;
	  }
	}
	
	// Add Sdigits
	if (GetDebug()>3) cerr<<"AliMUONSDigitzerv1::Exec TransientDigit to SDigit"<<endl;
	if ( digits[2] == icat ) muondata->AddSDigit(ich,tracks,charges,digits);
//	printf("test rm ich %d padX %d padY %d \n",ich, digits[0], digits[1]);
      }
      // Filling list of Sdigits per chamber for a given cathode.
      muondata->Fill("S");
      muondata->ResetSDigits();    
    } // end loop cathode
    fTDList->Delete();  
    
    for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
      if (fHitMap[ii]) {
	delete fHitMap[ii];
	fHitMap[ii] = 0;
      }
    }
    
    if (GetDebug()>2) 
      cerr<<"AliMUONSDigitizer::Exec: writing the TreeS: "
	  <<treeS->GetName()<<endl;

    runloader    = AliRunLoader::GetRunLoader();
    gime        = runloader->GetLoader("MUONLoader");
    gime->WriteSDigits("OVERWRITE");
    delete [] fHitMap;
    delete fTDList;
    muondata->ResetHits();
    gime->UnloadHits();
    gime->UnloadSDigits();
}
//------------------------------------------------------------------------
void AliMUONSDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
{
  //
  // Sort the list of tracks contributing to a given Sdigit
  // 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
d1775029	1
	2
	3	#include <Riostream.h>
	4	#include <TDirectory.h>
	5	#include <TFile.h>
	6	#include <TObjArray.h>
	7	#include <TPDGCode.h>
	8	#include <TTree.h>
	9	#include <TMath.h>
	10
	11	#include "AliMUON.h"
	12	#include "AliMUONChamber.h"
	13	#include "AliMUONConstants.h"
	14	#include "AliMUONDigit.h"
	15	#include "AliMUONSDigitizerv1.h"
	16	#include "AliMUONHit.h"
	17	#include "AliMUONHitMapA1.h"
	18	#include "AliMUONPadHit.h"
	19	#include "AliMUONTransientDigit.h"
	20	#include "AliRun.h"
	21	#include "AliRunLoader.h"
	22	#include "AliLoader.h"
	23
	24	ClassImp(AliMUONSDigitizerv1)
	25
	26	//___________________________________________
	27	AliMUONSDigitizerv1::AliMUONSDigitizerv1() :
	28	fHitMap(0),
	29	fTDList(0),
	30	fTDCounter(0),
	31	fDebug(0),
	32	fMask(0),
	33	fSignal(0)
	34	{
	35	// Default ctor - don't use it
	36	if (GetDebug()>2)
	37	cerr<<"AliMUONSDigitizerv1::AliMUONSDigitizerv1"
	38	<<"(AliRunDigitizer* manager) was processed"<<endl;
	39	}
	40
	41
	42	//------------------------------------------------------------------------
	43	AliMUONSDigitizerv1::~AliMUONSDigitizerv1()
	44	{
	45	// Destructor
	46	}
	47
	48	//------------------------------------------------------------------------
	49	void AliMUONSDigitizerv1::AddTransientDigit(AliMUONTransientDigit * mTD)
	50	{
	51	// Choosing the maping of the cathode plane of the chamber:
	52	Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
	53	fTDList->AddAtAndExpand(mTD, fTDCounter);
	54	fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
	55	fTDCounter++;
	56	}
	57
	58	//------------------------------------------------------------------------
	59	Bool_t AliMUONSDigitizerv1::ExistTransientDigit(AliMUONTransientDigit * mTD)
	60	{
	61	// Choosing the maping of the cathode plane of the chamber:
	62	Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
	63	return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
	64	}
65
66	//------------------------------------------------------------------------
67	Bool_t AliMUONSDigitizerv1::Init()
68	{
69
70	// Initialization
71	if (GetDebug()>2) Info("Init","AliMUONSDigitizerv1::Init() starts");
72
73	//Loaders (We assume input0 to be the output too)
74	AliRunLoader * runloader; // Input loader
75	AliLoader * gime;
76
77	// Getting runloader
78	runloader = AliRunLoader::GetRunLoader();
79	if (runloader == 0x0) {
80	Error("Init","RunLoader is not in input file 0");
81	return kFALSE; // RunDigitizer is not working.
82	}
83	// Getting MUONloader
84	gime = runloader->GetLoader("MUONLoader");
85	gime->LoadHits("READ");
86	gime->LoadSDigits("RECREATE");
87
88	return kTRUE;
89	}
90
91	//------------------------------------------------------------------------
92	void AliMUONSDigitizerv1::UpdateTransientDigit(Int_t track, AliMUONTransientDigit * mTD)
93	{
94	// Choosing the maping of the cathode plane of the chamber:
95	Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
96	AliMUONTransientDigit *pdigit =
97	static_cast<AliMUONTransientDigit*>(fHitMap[iNchCpl]->GetHit(mTD->PadX(),mTD->PadY()));
98	// update charge
99	//
100	Int_t iqpad = mTD->Signal(); // charge per pad
101	pdigit->AddSignal(iqpad);
102	pdigit->AddPhysicsSignal(iqpad);
103	// update list of tracks
104	//
105	Int_t charge;
106	track=+ fMask;
107	if (fSignal) charge = iqpad;
108	//else charge = kBgTag;
109	else charge = iqpad + fMask;
110
111	pdigit->UpdateTrackList(track,charge);
112	}
113
114
115	//--------------------------------------------------------------------------
116	void AliMUONSDigitizerv1::MakeTransientDigit(Int_t track, Int_t iHit, AliMUONHit * mHit)
117	{
118	AliMUON pMUON = (AliMUON ) gAlice->GetModule("MUON");
119	if (!pMUON) {
120	cerr<<"AliMUONSDigitizerv1::MakeTransientDigit Error:"
121	<<" module MUON not found in the input file"<<endl;
122	}
123	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit starts"<<endl;
124	Int_t ichamber = mHit->Chamber()-1;
125	AliMUONChamber & chamber = pMUON->Chamber(ichamber);
126	Float_t xhit = mHit->X();
127	Float_t yhit = mHit->Y();
128	Float_t zhit = mHit->Z();
129	Float_t eloss= mHit->Eloss();
130	Float_t tof = mHit->Age();
131	// Variables for chamber response from AliMUONChamber::DisIntegration
132	Float_t newdigit[6][500]; // Pad information
133	Int_t nnew=0; // Number of touched Pads per hit
134	Int_t digits[6];
135
136	//
137	// Calls the charge disintegration method of the current chamber
138	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit calling AliMUONChamber::DisIngtegration starts"<<endl;
139	chamber.DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newdigit);
140	// Creating a new TransientDigits from hit
141	for(Int_t iTD=0; iTD<nnew; iTD++) {
142	digits[0] = Int_t(newdigit[1][iTD]); // Padx of the Digit
143	digits[1] = Int_t(newdigit[2][iTD]); // Pady of the Digit
144	digits[2] = Int_t(newdigit[5][iTD]); // Cathode plane
145	digits[3] = Int_t(newdigit[3][iTD]); // Induced charge in the Pad !!!! Int_t not correct
146	if (fSignal) digits[4] = Int_t(newdigit[3][iTD]);
147	else digits[4] = 0;
148	digits[5] = iHit+fMask; // Hit number in the list
149	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit " <<
150	cerr<<"AliMUONSDigitizerv1::MakeTransientDigit " <<
151	"PadX "<< digits[0] << " " <<
152	"PadY "<< digits[1] << " " <<
153	"Plane " << digits[2] << " " <<
154	"Charge " << digits[3] <<" " <<
155	"newdigit " << newdigit[3][iTD] <<" " <<
156	"Hit " << digits[5] << endl;
157	// list of tracks
158	Int_t charge;
159	track += fMask;
160	if (fSignal) charge = digits[3];
161	//else charge = kBgTag;
162	else charge = digits[3] + fMask;
163
164	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit Creating AliMUONTransientDigit"<<endl;
165	AliMUONTransientDigit * mTD = new AliMUONTransientDigit(ichamber, digits);
166	mTD->AddToTrackList(track,charge);
167	if (!ExistTransientDigit(mTD)) {
168	AddTransientDigit(mTD);
169	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit Adding TransientDigit"<<endl;
170	}
171	else {
172	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::MakeTransientDigit updating TransientDigit"<<endl;
173	UpdateTransientDigit(track, mTD);
174	delete mTD;
175	}
176	}
177	}
178	//-----------------------------------------------------------------------
179	void AliMUONSDigitizerv1::Exec(Option_t* option)
180	{
181	TString optionString = option;
182	if (optionString.Data() == "deb") {
183	Info("SDigitize","Called with option deb ");
184	fDebug = 3;
185	}
186
187	AliMUONChamber* chamber;
188	AliSegmentation* c1Segmentation; //Cathode plane c1 of the chamber
189	AliSegmentation* c2Segmentation; //Cathode place c2 of the chamber
190
191	if (GetDebug()>2) Info("SDigitize","AliMUONSDigitizerv1::Exec starts");
192	fTDList = new TObjArray;
193
194	//Loaders (We assume input0 to be the output too)
195	AliRunLoader * runloader; // Input loader
196	AliLoader * gime;
197
198	// Getting runloader
199	runloader = AliRunLoader::GetRunLoader();
200	if (runloader == 0x0) {
201	Error("SDigitize","RunLoader is not in input file 0");
202	return; // RunDigitizer is not working.
203	}
204	// Getting MUONloader
205	gime = runloader->GetLoader("MUONLoader");
206	if (gime->TreeH()==0x0) {
207	if (GetDebug()>2) Info("SDigitize","TreeH is not loaded yet. Loading...");
208	gime->LoadHits("READ");
209	if (GetDebug()>2) Info("SDigitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
210	}
211
212	if (GetDebug()>2) Info("SDigitize","Loaders ready");
213
214	if (runloader->GetAliRun() == 0x0) runloader->LoadgAlice();
215	gAlice = runloader->GetAliRun();
216
217	// Getting Module MUON
218	AliMUON pMUON = (AliMUON ) gAlice->GetDetector("MUON");
219	if (!pMUON) {
220	Error("SDigitize","Module MUON not found in the input file");
221	return;
222	}
223	// Getting Muon data
224	AliMUONData * muondata = pMUON->GetMUONData();
225	muondata->SetLoader(gime);
226	muondata->SetTreeAddress("H");
227
228	Int_t currentevent = runloader->GetEventNumber();
229
230	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::Exec Event Number is "<<currentevent <<endl;
231	if ( (currentevent<10) \|\|
232	(Int_t(TMath::Log10(currentevent)) == TMath::Log10(currentevent) ) )
233	cout <<"AliMUONSDigitizerv1::Exec Event Number is "<< currentevent <<endl;
234
235	// New branch per chamber for MUON digit in the tree of digits
236	if (gime->TreeS() == 0x0) {
237	gime->MakeSDigitsContainer();
238	}
239	TTree* treeS = gime->TreeS();
240	muondata->MakeBranch("S");
241	muondata->SetTreeAddress("S");
242
243	// Array of pointer of the AliMUONHitMapA1:
244	// two HitMaps per chamber, or one HitMap per cahtode plane
245	fHitMap= new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];
246
247	//Loop over chambers for the definition AliMUONHitMap
248	for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
249	chamber = &(pMUON->Chamber(i));
250	c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
251	fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
252	c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
253	fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
254	}
255
256	// Loop to Sdigitize
257	fSignal = kTRUE;
258
259	// Setting the address
260	TTree *treeH = gime->TreeH();
261	if (treeH == 0x0) {
262	Error("SDigitize","Can not get TreeH from input ");
263	Info("SDigitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
264	return;
265	}
266	if (GetDebug()>2) {
267	cerr<<"AliMUONSDigitizerv1::Exec treeH" << treeH <<endl;
268	}
269
270	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::Exec Setting tree addresses"<<endl;
271
272	//
273	// Loop over tracks
274	Int_t itrack;
275	Int_t ntracks = (Int_t) treeH->GetEntries();
276	for (itrack = 0; itrack < ntracks; itrack++) {
277	if (GetDebug()>2) cerr<<"AliMUONSDigitizerv1::Exec itrack = "<<itrack<<endl;
278	muondata->ResetHits();
279	treeH->GetEvent(itrack);
280	//
281	// Loop over hits
282	Int_t ihit, ichamber;
283	AliMUONHit* mHit;
284	TClonesArray* hits = muondata->Hits();
285	for(ihit = 0; ihit < hits->GetEntriesFast(); ihit++) {
286	mHit = static_cast<AliMUONHit*>(hits->At(ihit));
287	ichamber = mHit->Chamber()-1; // chamber number
288	if (ichamber > AliMUONConstants::NCh()-1) {
289	cerr<<"AliMUONSDigitizer: ERROR: "
290	<<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
291	<<ichamber<<", "<< AliMUONConstants::NCh()<<endl;
292	return;
293	}
294	chamber = &(pMUON->Chamber(ichamber));
295	//
296	//Dumping Hit content:
297	if (GetDebug()>2) {
298	cerr<<"AliMuonDigitizerv1::Exec ihit, ichamber, x, y, z, eloss " <<
299	ihit << " " <<
300	mHit->Chamber() << " " <<
301	mHit->X() << " " <<
302	mHit->Y() << " " <<
303	mHit->Z() << " " <<
304	mHit->Eloss() << " " << endl;
305	}
306	//
307	// Inititializing Correlation
308	chamber->ChargeCorrelationInit();
309	if (ichamber < AliMUONConstants::NTrackingCh()) {
310	// Tracking Chamber
311	// Initialize hit position (cursor) in the segmentation model
312	chamber->SigGenInit(mHit->X(), mHit->Y(), mHit->Z());
313	} else {
314	// Trigger Chamber
315	}
316	MakeTransientDigit(itrack, ihit, mHit);
317	} // hit loop
318	} // track loop
319	if (GetDebug()>2) cerr<<"AliMUONSDigitizer::Exec End of hits, track and file loops"<<endl;
320
321	// Loop on cathodes
322	Int_t icat;
323	for(icat=0; icat<2; icat++) {
324	//
325	// Filling SDigit List
326	Int_t tracks[kMAXTRACKS];
327	Int_t charges[kMAXTRACKS];
328	Int_t nentries = fTDList->GetEntriesFast();
329	Int_t digits[6];
330	for (Int_t nent = 0; nent < nentries; nent++) {
331	AliMUONTransientDigit address = (AliMUONTransientDigit)fTDList->At(nent);
332	if (address == 0) continue;
333	Int_t ich = address->Chamber();
334	Int_t q = address->Signal();
335
336	if (!q) continue; // not mandatory ?
337
338	digits[0] = address->PadX();
339	digits[1] = address->PadY();
340	digits[2] = address->Cathode()-1;
341	digits[3] = q;
342	digits[4] = address->Physics();
343	digits[5] = address->Hit();
344
345	Int_t nptracks = address->GetNTracks();
346
347	if (nptracks > kMAXTRACKS) {
348	if (GetDebug() >0) {
349	cerr<<"AliMUONSDigitizer:Exec nptracks > 10 "<<nptracks;
350	cerr<<"reset to max value "<<kMAXTRACKS<<endl;
351	}
352	nptracks = kMAXTRACKS;
353	}
354	if (nptracks > 2 && GetDebug() >2) {
355	cerr<<"AliMUONSDigitizer::Exec nptracks > 2 "<<nptracks<<endl;
356	// printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
357	}
358	for (Int_t tr = 0; tr < nptracks; tr++) {
359	tracks[tr] = address->GetTrack(tr);
360	charges[tr] = address->GetCharge(tr);
361	} //end loop over list of tracks for one pad
362	// Sort list of tracks according to charge
363	if (nptracks > 1) {
364	SortTracks(tracks,charges,nptracks);
365	}
366	if (nptracks < kMAXTRACKS ) {
367	for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
368	tracks[i] = -1;
369	charges[i] = 0;
370	}
371	}
372
373	// Add Sdigits
374	if (GetDebug()>3) cerr<<"AliMUONSDigitzerv1::Exec TransientDigit to SDigit"<<endl;
375	if ( digits[2] == icat ) muondata->AddSDigit(ich,tracks,charges,digits);
376	// printf("test rm ich %d padX %d padY %d \n",ich, digits[0], digits[1]);
377	}
378	// Filling list of Sdigits per chamber for a given cathode.
379	muondata->Fill("S");
380	muondata->ResetSDigits();
381	} // end loop cathode
382	fTDList->Delete();
383
384	for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
385	if (fHitMap[ii]) {
386	delete fHitMap[ii];
387	fHitMap[ii] = 0;
388	}
389	}
390
391	if (GetDebug()>2)
392	cerr<<"AliMUONSDigitizer::Exec: writing the TreeS: "
393	<<treeS->GetName()<<endl;
394
395	runloader = AliRunLoader::GetRunLoader();
396	gime = runloader->GetLoader("MUONLoader");
397	gime->WriteSDigits("OVERWRITE");
398	delete [] fHitMap;
399	delete fTDList;
400	muondata->ResetHits();
401	gime->UnloadHits();
402	gime->UnloadSDigits();
403	}
404	//------------------------------------------------------------------------
405	void AliMUONSDigitizerv1::SortTracks(Int_t tracks,Int_t charges,Int_t ntr)
406	{
407	//
408	// Sort the list of tracks contributing to a given Sdigit
409	// Only the 3 most significant tracks are acctually sorted
410	//
411
412	//
413	// Loop over signals, only 3 times
414	//
415
416	Int_t qmax;
417	Int_t jmax;
418	Int_t idx[3] = {-2,-2,-2};
419	Int_t jch[3] = {-2,-2,-2};
420	Int_t jtr[3] = {-2,-2,-2};
421	Int_t i,j,imax;
422
423	if (ntr<3) imax=ntr;
424	else imax=3;
425	for(i=0;i<imax;i++){
426	qmax=0;
427	jmax=0;
428
429	for(j=0;j<ntr;j++){
430
431	if((i == 1 && j == idx[i-1])
432	\|\|(i == 2 && (j == idx[i-1] \|\| j == idx[i-2]))) continue;
433
434	if(charges[j] > qmax) {
435	qmax = charges[j];
436	jmax=j;
437	}
438	}
439
440	if(qmax > 0) {
441	idx[i]=jmax;
442	jch[i]=charges[jmax];
443	jtr[i]=tracks[jmax];
444	}
445
446	}
447
448	for(i=0;i<3;i++){
449	if (jtr[i] == -2) {
450	charges[i]=0;
451	tracks[i]=0;
452	} else {
453	charges[i]=jch[i];
454	tracks[i]=jtr[i];
455	}
456	}
457	}