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


//  Do the Digitization (Digit) from summable Digits (SDigit)
//  Allow the merging of signal file with background file(s).

#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 "AliMUONDigitizerv2.h"
#include "AliMUONHit.h"
#include "AliMUONHitMapA1.h"
#include "AliMUONPadHit.h"
#include "AliMUONTransientDigit.h"
#include "AliRun.h"
#include "AliRunDigitizer.h"
#include "AliRunLoader.h"
#include "AliLoader.h"

ClassImp(AliMUONDigitizerv2)

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

//___________________________________________
AliMUONDigitizerv2::AliMUONDigitizerv2(AliRunDigitizer* manager):
  AliDigitizer(manager),
  fHitMap(0),
  fTDList(0),
  fTDCounter(0),
  fDebug(0),
  fMask(0),
  fSignal(0)
{
// ctor which should be used
  if (GetDebug()>2) 
    cerr<<"AliMUONDigitizerv2::AliMUONDigitizerv2"
	<<"(AliRunDigitizer* manager) was processed"<<endl;
}

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

//------------------------------------------------------------------------
void AliMUONDigitizerv2::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 AliMUONDigitizerv2::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 AliMUONDigitizerv2::Init()
{

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

  //Loaders (We assume input0 to be the output too)
  AliRunLoader * runloader;    // Input   loader
  AliLoader    * gime; 
 
  // Getting runloader
  runloader    = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
  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->LoadSDigits("READ");
  gime->LoadDigits("RECREATE");
  
  return kTRUE;
}

//------------------------------------------------------------------------
void AliMUONDigitizerv2::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
  //
  // List of tracks and trackcharge
  Int_t itrack;
  for(itrack=0;itrack<kMAXTRACKS;itrack++) {
    pdigit->UpdateTrackList(mTD->Track(itrack), mTD->TrackCharge(itrack));
  }
}

//--------------------------------------------------------------------------
void AliMUONDigitizerv2::MakeTransientDigitFromSDigit(Int_t iChamber, AliMUONDigit * sDigit)
{
  Int_t digits[6];
  

  //
  // Creating a new TransientDigits from SDigit
  digits[0] = sDigit->PadX();  // Padx of the Digit
  digits[1] = sDigit->PadY();  // Pady of the Digit
  digits[2] = sDigit->Cathode()+1;  // Cathode plane
  digits[3] = sDigit->Signal();  // Induced charge in the Pad
  if (fSignal) digits[4] = sDigit->Signal();
  else         digits[4] = 0;
  digits[5] = sDigit->Hit();    // Hit number in the list
  if (GetDebug()>2) cerr<<"AliMUONDigitizerv2::MakeTransientDigitFromSDigit " <<
		      "PadX "<< digits[0] << " " <<
		      "PadY "<< digits[1] << " " <<
		      "Plane " << digits[2] << " " <<
		      "Charge " << digits[3] <<" " <<
		      "Hit " << digits[5] << endl;   
		 
  if (GetDebug()>2) cerr<<"AliMUONDigitizerv2::MakeTransientDigitFromSDigit Creating AliMUONTransientDigit"<<endl;
  AliMUONTransientDigit * mTD = new AliMUONTransientDigit(iChamber, digits);

  // List of tracks and trackcharge
  Int_t itrack;
  for(itrack=0;itrack<kMAXTRACKS;itrack++) {
    mTD->AddToTrackList(fMask+sDigit->Track(itrack), sDigit->TrackCharge(itrack));
  }

  if (!ExistTransientDigit(mTD)) {
    AddTransientDigit(mTD);
    if (GetDebug()>2) cerr<<"AliMUONDigitizerv2::MakeTransientDigitFromSDigit Adding TransientDigit"<<endl;
  }
  else {
    if (GetDebug()>2) cerr<<"AliMUONDigitizerv2::MakeTransientDigitFromSDigit updating TransientDigit"<<endl;
    UpdateTransientDigit(0, mTD);
    delete mTD;
  }
}

//-----------------------------------------------------------------------
void AliMUONDigitizerv2::Exec(Option_t* option)
{
  TString optionString = option;
  if (optionString.Data() == "deb") {
    Info("Digitize","Called with option deb ");
    fDebug = 3;
  }

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

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

  //Loaders (We assume input0 to be the output too)
  AliRunLoader *runloader, *runloaderOut;    // Input /  Output loaders
  AliLoader    *gime, *gimeOut; 
 
  // Getting runloader
  runloader    = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
  if (runloader == 0x0) {
    Error("Exec","RunLoader is not in input file 0");
    return;  // RunDigitizer is not working.  
  }
  // Getting MUONloader
  gime        = runloader->GetLoader("MUONLoader");
  if (gime->TreeS()==0x0) {
      if (GetDebug()>2) Info("Exec","TreeS is not loaded yet. Loading...");
     gime->LoadSDigits("READ");
       if (GetDebug()>2) Info("Exec","Now treeS is %#x. MUONLoader is %#x",gime->TreeS(),gime);
  }

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

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

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

  // Loading Event
  Int_t currentevent = fManager->GetOutputEventNr();
  
  if (GetDebug()>2) cerr<<"AliMUONDigitizerv2::Exec() Event Number is "<<currentevent <<endl;
  if ( (currentevent<10)                                                 || 
       (Int_t(TMath::Log10(currentevent)) == TMath::Log10(currentevent) ) )
    cout <<"ALiMUONDigitizerv2::Exec() Event Number is "<< currentevent <<endl;

  // Output runloader 
  runloaderOut  = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
  gimeOut = runloaderOut->GetLoader("MUONLoader");
  

  if (gimeOut->TreeD() == 0x0) {
    gimeOut->MakeDigitsContainer();
  }
  TTree* treeD = gimeOut->TreeD();
  muondata->SetLoader(gimeOut);
  muondata->MakeBranch("D");
  muondata->SetTreeAddress("D");

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