[u/mrichter/AliRoot.git] / RICH / AliRICHMerger.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/10/14 14:57:32  hristov
Merging the VirtualMC branch to the main development branch (HEAD)

Revision 1.7.6.1  2002/07/24 10:07:52  alibrary
Updating VirtualMC

Revision 1.7  2001/11/02 15:37:26  hristov
Digitizer class created. Code cleaning and bug fixes (J.Chudoba)

Revision 1.5  2001/10/23 13:03:35  hristov
The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba)

Revision 1.4  2001/10/21 18:31:24  hristov
Several pointers were set to zero in the default constructors to avoid memory management problems

Revision 1.3  2001/05/16 14:57:20  alibrary
New files for folders and Stack

Revision 1.2  2001/03/14 18:16:08  jbarbosa
Corrected bug (more to correct).
File "points.dat" is no longer created.

Revision 1.1  2001/02/27 22:13:34  jbarbosa
Implementing merger class.

*/
#include <Riostream.h> 

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

#include "AliRICHMerger.h"
#include "AliRICHChamber.h"
#include "AliHitMap.h"
#include "AliRICHHitMapA1.h"
#include "AliRICH.h"
#include "AliRICHHit.h"
#include "AliRICHSDigit.h"
#include "AliRICHDigit.h"
#include "AliRICHTransientDigit.h"
#include "AliRun.h"
#include "AliPDG.h"

ClassImp(AliRICHMerger)

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

//------------------------------------------------------------------------
AliRICHMerger::~AliRICHMerger()
{
// Destructor
  if (fTrH1)       delete fTrH1;
  if (fHitsBgr)    delete fHitsBgr;
  if (fSDigitsBgr) delete fSDigitsBgr;
  if (fHitMap)     delete fHitMap;
  if (fList)       delete fList;
  if (fBgrFile)    delete fBgrFile;
}

//------------------------------------------------------------------------
Bool_t AliRICHMerger::Exists(const AliRICHSDigit *padhit)
{
  return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
}

//------------------------------------------------------------------------
void AliRICHMerger::Update(AliRICHSDigit *padhit)
{
  AliRICHTransientDigit *pdigit = 
    static_cast<AliRICHTransientDigit*>(
      fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));

  // update charge
  //
  Int_t iqpad    = Int_t(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 AliRICHMerger::CreateNew(AliRICHSDigit *padhit)
{
  fList->AddAtAndExpand(
    new AliRICHTransientDigit(fNch,fDigits),fCounter);
  fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);

  AliRICHTransientDigit* pdigit = 
    static_cast<AliRICHTransientDigit*>(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 AliRICHMerger::Init()
{
// Initialisation
    
  if (fMerge && !fBgrFile) fBgrFile = InitBgr();
}


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

//------------------------------------------------------------------------
void AliRICHMerger::Digitise(Int_t nev, Int_t flag)
{

  // keep galice.root for signal and name differently the file for 
  // background when add! otherwise the track info for signal will be lost !

  Int_t particle;

  AliRICHChamber*       iChamber;
  AliSegmentation*  segmentation;
  
  fList = new TObjArray;
    
  AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
  fHitMap= new AliHitMap* [kNCH];
    
  if (fMerge ) {
    fBgrFile->cd();
    //
    // Get Hits Tree header from file
    if(fHitsBgr) fHitsBgr->Clear();
    if(fSDigitsBgr) fSDigitsBgr->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",pRICH->GetName());
    if (fTrH1 && fHitsBgr) {
      branch = fTrH1->GetBranch(branchname);
      if (branch) branch->SetAddress(&fHitsBgr);
    }
    if (fTrH1 && fSDigitsBgr) {
      branch = fTrH1->GetBranch("RICHSDigits");
      if (branch) branch->SetAddress(&fSDigitsBgr);
    }
  }
    
  for (Int_t i =0; i<kNCH; i++) {
    iChamber= &(pRICH->Chamber(i));
    segmentation=iChamber->GetSegmentationModel(1);
    fHitMap[i] = new AliRICHHitMapA1(segmentation, fList);
  }
  //
  //   Loop over tracks
  //
    
  fSignal = kTRUE;
  fCounter = 0;
  TTree *treeH = gAlice->TreeH();
  Int_t ntracks =(Int_t) treeH->GetEntries();
  for (fTrack=0; fTrack<ntracks; fTrack++) {
    gAlice->ResetHits();
    treeH->GetEvent(fTrack);
    //
    //   Loop over hits
    for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1); 
	mHit;
	mHit=(AliRICHHit*)pRICH->NextHit()) 
    {
	  
      fNch = mHit->Chamber()-1;  // chamber number
      Int_t trackIndex = mHit->Track();
      if (fNch >= kNCH) {
	cerr<<"AliRICHMerger: chamber nr. fNch out of range: "<<fNch<<endl;
	cerr<<"               track: "<<fTrack<<endl;
	continue;
      }
      iChamber = &(pRICH->Chamber(fNch));
	  

// 
// If flag is set, create digits only for some particles
//
      TParticle *current = (TParticle*)gAlice->Particle(trackIndex);
      if (current->GetPdgCode() >= 50000050)
      {
	TParticle *motherofcurrent = (TParticle*)gAlice->Particle(current->GetFirstMother());
	particle = motherofcurrent->GetPdgCode();
      }
      else
      {
	particle = current->GetPdgCode();
      }

	  
      //printf("Flag:%d\n",flag);
      //printf("Track:%d\n",track);
      //printf("Particle:%d\n",particle);
	  
      Int_t digitise=1;
	  
      if (flag == 1) 
	if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
	  digitise=0;
	  
      if (flag == 2)
	if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
	   || TMath::Abs(particle)==311)
	  digitise=0;
	  
      if (flag == 3 && TMath::Abs(particle)==2212)
	digitise=0;
	  
      if (flag == 4 && TMath::Abs(particle)==13)
	digitise=0;
	  
      if (flag == 5 && TMath::Abs(particle)==11)
	digitise=0;
	  
      if (flag == 6 && TMath::Abs(particle)==2112)
	digitise=0;
	  
	  
      //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise); 
	  
	  
      if (digitise)
      {
	      
	//
	// Loop over pad hits
	for (AliRICHSDigit* mPad=
	       (AliRICHSDigit*)pRICH->FirstPad(mHit,pRICH->SDigits());
	     mPad;
	     mPad=(AliRICHSDigit*)pRICH->NextPad(pRICH->SDigits()))
	{
	  Int_t iqpad    = mPad->QPad();       // charge per pad
	  fDigits[0]=mPad->PadX();
	  fDigits[1]=mPad->PadY();
	  fDigits[2]=iqpad;
	  fDigits[3]=iqpad;
	  fDigits[4]=mPad->HitNumber();
		  
	  // build the list of fired pads and update the info
	  if (Exists(mPad)) {
	    Update(mPad);
	  } else {
	    CreateNew(mPad);
	  }
	} //end loop over clusters
      }// track type condition
    } // 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 (fSDigitsBgr)    fSDigitsBgr->Clear();
	
      fTrH1->GetEvent(fTrack);
      //
      //   Loop over hits
      AliRICHHit* mHit;
      for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i) 
      {	
	mHit   = (AliRICHHit*) (*fHitsBgr)[i];
	fNch   = mHit->Chamber()-1;  // chamber number
	iChamber = &(pRICH->Chamber(fNch));

	//
	// Loop over pad hits
	for (AliRICHSDigit* mPad =
	       (AliRICHSDigit*)pRICH->FirstPad(mHit,fSDigitsBgr);
	     mPad;
	     mPad = (AliRICHSDigit*)pRICH->NextPad(fSDigitsBgr))
	{
	  fDigits[0] = mPad->PadX(); 
	  fDigits[1] = mPad->PadY();  
	  fDigits[2] = mPad->QPad();
	  fDigits[3] = 0;
	  fDigits[4] = -1; // set hit number to -1 for bgr
		
	  // 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[kMAXTRACKSPERRICHDIGIT];
  Int_t charges[kMAXTRACKSPERRICHDIGIT];
  Int_t nentries=fList->GetEntriesFast();
    
  for (Int_t nent=0;nent<nentries;nent++) {
    AliRICHTransientDigit *transDigit=(AliRICHTransientDigit*)fList->At(nent);
    if (transDigit==0) continue; 
    Int_t ich=transDigit->GetChamber();
    Int_t q=transDigit->Signal(); 
    iChamber=&(pRICH->Chamber(ich));
    AliRICHResponse * response=iChamber->GetResponseModel();
    Int_t adcmax= (Int_t) response->MaxAdc();
      
      
    // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
    //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
    Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());

    //printf("Pedestal:%d\n",pedestal);
    //Int_t pedestal=0;
    Float_t treshold = (pedestal + 4*2.2);
      
    Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
    Float_t noise     = gRandom->Gaus(0, meanNoise);
    q+=(Int_t)(noise + pedestal);
    //q+=(Int_t)(noise);
    //          magic number to be parametrised !!! 
    if ( q <= treshold) 
    {
      q = q - pedestal;
      continue;
    }
    q = q - pedestal;
    if ( q >= adcmax) q=adcmax;
    fDigits[0]=transDigit->PadX();
    fDigits[1]=transDigit->PadY();
    fDigits[2]=q;
    fDigits[3]=transDigit->Physics();
    fDigits[4]=transDigit->Hit();

    Int_t nptracks = transDigit->GetNTracks();  

    // this was changed to accomodate the real number of tracks
    if (nptracks > kMAXTRACKSPERRICHDIGIT) {
      printf("Attention - tracks > 10 %d\n",nptracks);
      nptracks=kMAXTRACKSPERRICHDIGIT;
    }
    if (nptracks > 2) {
      printf("Attention - tracks > 2  %d \n",nptracks);
      //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]=transDigit->GetTrack(tr);
      charges[tr]=transDigit->GetCharge(tr);
    }      //end loop over list of tracks for one pad
    if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
      for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
	tracks[t]=0;
	charges[t]=0;
      }
    }
    //write file
    //if (ich==2)
    //fprintf(points,"%4d,      %4d,      %4d\n",digits[0],digits[1],digits[2]);
      
    // fill digits
    pRICH->AddDigits(ich,tracks,charges,fDigits);
  }	
  gAlice->TreeD()->Fill();

  fList->Delete();
  for(Int_t ii=0;ii<kNCH;++ii) {
    if (fHitMap[ii]) {
      delete fHitMap[ii];
      fHitMap[ii]=0;
    }
  }
    
  TClonesArray *richDigits;
  for (Int_t k=0;k<kNCH;k++) {
    richDigits = pRICH->DigitsAddress(k);
    Int_t ndigit=richDigits->GetEntriesFast();
    printf ("Chamber %d digits %d \n",k,ndigit);
  }
  pRICH->ResetDigits(); /// ??? should it be here???
  gAlice->TreeD()->Write(0,TObject::kOverwrite);
  // reset tree
  //    gAlice->TreeD()->Reset();
  //    delete fList; // deleted in dtor
  // gObjectTable->Print();

}
Commit	Line	Data
d8a72780	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$
7118aef0	18	Revision 1.9 2002/10/14 14:57:32 hristov
	19	Merging the VirtualMC branch to the main development branch (HEAD)
	20
b9d0a01d	21	Revision 1.7.6.1 2002/07/24 10:07:52 alibrary
	22	Updating VirtualMC
	23
adf7bb67	24	Revision 1.7 2001/11/02 15:37:26 hristov
	25	Digitizer class created. Code cleaning and bug fixes (J.Chudoba)
	26
31eeed4f	27	Revision 1.5 2001/10/23 13:03:35 hristov
	28	The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba)
	29
6e585aa2	30	Revision 1.4 2001/10/21 18:31:24 hristov
	31	Several pointers were set to zero in the default constructors to avoid memory management problems
	32
2685bf00	33	Revision 1.3 2001/05/16 14:57:20 alibrary
	34	New files for folders and Stack
	35
9e1a0ddb	36	Revision 1.2 2001/03/14 18:16:08 jbarbosa
	37	Corrected bug (more to correct).
	38	File "points.dat" is no longer created.
	39
14886cb7	40	Revision 1.1 2001/02/27 22:13:34 jbarbosa
	41	Implementing merger class.
	42
d8a72780	43	*/
7118aef0	44	#include <Riostream.h>
d8a72780	45
d8a72780	46	#include <TTree.h>
d8a72780	47	#include <TObjArray.h>
	48	#include <TFile.h>
	49	#include <TDirectory.h>
	50	#include <TParticle.h>
	51
d8a72780	52	#include "AliRICHMerger.h"
	53	#include "AliRICHChamber.h"
	54	#include "AliHitMap.h"
	55	#include "AliRICHHitMapA1.h"
	56	#include "AliRICH.h"
	57	#include "AliRICHHit.h"
	58	#include "AliRICHSDigit.h"
	59	#include "AliRICHDigit.h"
	60	#include "AliRICHTransientDigit.h"
	61	#include "AliRun.h"
	62	#include "AliPDG.h"
	63
	64	ClassImp(AliRICHMerger)
	65
	66	//___________________________________________
b762c2f6	67	AliRICHMerger::AliRICHMerger()
d8a72780	68	{
d8a72780	69	// Default constructor
b762c2f6	70	fEvNrSig = 0;
	71	fEvNrBgr = 0;
	72	fMerge = kDigitize;
	73	fFnBgr = 0;
	74	fHitMap = 0;
	75	fList = 0;
	76	fTrH1 = 0;
	77	fHitsBgr = 0;
	78	fSDigitsBgr = 0;
	79	fHitMap = 0;
	80	fList = 0;
	81	fBgrFile = 0;
d8a72780	82	}
	83
	84	//------------------------------------------------------------------------
	85	AliRICHMerger::~AliRICHMerger()
	86	{
	87	// Destructor
b762c2f6	88	if (fTrH1) delete fTrH1;
	89	if (fHitsBgr) delete fHitsBgr;
	90	if (fSDigitsBgr) delete fSDigitsBgr;
	91	if (fHitMap) delete fHitMap;
	92	if (fList) delete fList;
	93	if (fBgrFile) delete fBgrFile;
d8a72780	94	}
	95
	96	//------------------------------------------------------------------------
b762c2f6	97	Bool_t AliRICHMerger::Exists(const AliRICHSDigit *padhit)
d8a72780	98	{
b762c2f6	99	return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
d8a72780	100	}
	101
	102	//------------------------------------------------------------------------
b762c2f6	103	void AliRICHMerger::Update(AliRICHSDigit *padhit)
d8a72780	104	{
b762c2f6	105	AliRICHTransientDigit *pdigit =
	106	static_cast<AliRICHTransientDigit*>(
	107	fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
	108
	109	// update charge
	110	//
	111	Int_t iqpad = Int_t(padhit->QPad()); // charge per pad
	112	pdigit->AddSignal(iqpad);
	113	pdigit->AddPhysicsSignal(iqpad);
	114
	115	// update list of tracks
	116	//
	117	Int_t track, charge;
	118	if (fSignal) {
	119	track = fTrack;
	120	charge = iqpad;
	121	} else {
	122	track = kBgTag;
	123	charge = kBgTag;
	124	}
	125	pdigit->UpdateTrackList(track,charge);
d8a72780	126	}
	127
	128	//------------------------------------------------------------------------
b762c2f6	129	void AliRICHMerger::CreateNew(AliRICHSDigit *padhit)
d8a72780	130	{
b762c2f6	131	fList->AddAtAndExpand(
	132	new AliRICHTransientDigit(fNch,fDigits),fCounter);
	133	fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
	134
	135	AliRICHTransientDigit* pdigit =
	136	static_cast<AliRICHTransientDigit*>(fList->Last());
	137	// list of tracks
	138	Int_t track, charge;
	139	if (fSignal) {
	140	track = fTrack;
	141	charge = padhit->QPad();
	142	} else {
	143	track = kBgTag;
	144	charge = kBgTag;
	145	}
	146	pdigit->AddToTrackList(track,charge);
	147	fCounter++;
d8a72780	148	}
	149
	150
	151	//------------------------------------------------------------------------
	152	void AliRICHMerger::Init()
	153	{
	154	// Initialisation
	155
b762c2f6	156	if (fMerge && !fBgrFile) fBgrFile = InitBgr();
d8a72780	157	}
	158
	159
	160
	161	//------------------------------------------------------------------------
	162	TFile* AliRICHMerger::InitBgr()
	163	{
	164	// Initialise background event
b762c2f6	165	TFile *file = new TFile(fFnBgr);
d8a72780	166	// add error checking later
b762c2f6	167	printf("\n AliRICHMerger has opened %s file with background event \n", fFnBgr);
	168	fHitsBgr = new TClonesArray("AliRICHHit",1000);
	169	fSDigitsBgr = new TClonesArray("AliRICHSDigit",1000);
	170	return file;
d8a72780	171	}
	172
	173	//------------------------------------------------------------------------
	174	void AliRICHMerger::Digitise(Int_t nev, Int_t flag)
	175	{
	176
b762c2f6	177	// keep galice.root for signal and name differently the file for
b762c2f6	178	// background when add! otherwise the track info for signal will be lost !
d8a72780	179
	180	Int_t particle;
	181
d8a72780	182	AliRICHChamber* iChamber;
	183	AliSegmentation* segmentation;
	184
b762c2f6	185	fList = new TObjArray;
d8a72780	186
b762c2f6	187	AliRICH pRICH = (AliRICH ) gAlice->GetDetector("RICH");
b762c2f6	188	fHitMap= new AliHitMap* [kNCH];
d8a72780	189
b762c2f6	190	if (fMerge ) {
	191	fBgrFile->cd();
	192	//
	193	// Get Hits Tree header from file
	194	if(fHitsBgr) fHitsBgr->Clear();
	195	if(fSDigitsBgr) fSDigitsBgr->Clear();
	196	if(fTrH1) delete fTrH1;
	197	fTrH1 = 0;
d8a72780	198
b762c2f6	199	char treeName[20];
	200	sprintf(treeName,"TreeH%d",fEvNrBgr);
	201	fTrH1 = (TTree*)gDirectory->Get(treeName);
	202	if (!fTrH1) {
	203	printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
d8a72780	204	}
d8a72780	205	//
b762c2f6	206	// Set branch addresses
	207	TBranch *branch;
	208	char branchname[20];
	209	sprintf(branchname,"%s",pRICH->GetName());
	210	if (fTrH1 && fHitsBgr) {
	211	branch = fTrH1->GetBranch(branchname);
	212	if (branch) branch->SetAddress(&fHitsBgr);
	213	}
	214	if (fTrH1 && fSDigitsBgr) {
	215	branch = fTrH1->GetBranch("RICHSDigits");
	216	if (branch) branch->SetAddress(&fSDigitsBgr);
	217	}
	218	}
d8a72780	219
b762c2f6	220	for (Int_t i =0; i<kNCH; i++) {
	221	iChamber= &(pRICH->Chamber(i));
	222	segmentation=iChamber->GetSegmentationModel(1);
	223	fHitMap[i] = new AliRICHHitMapA1(segmentation, fList);
	224	}
	225	//
	226	// Loop over tracks
	227	//
	228
	229	fSignal = kTRUE;
	230	fCounter = 0;
	231	TTree *treeH = gAlice->TreeH();
	232	Int_t ntracks =(Int_t) treeH->GetEntries();
	233	for (fTrack=0; fTrack<ntracks; fTrack++) {
	234	gAlice->ResetHits();
	235	treeH->GetEvent(fTrack);
	236	//
	237	// Loop over hits
	238	for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
	239	mHit;
	240	mHit=(AliRICHHit*)pRICH->NextHit())
	241	{
d8a72780	242
b762c2f6	243	fNch = mHit->Chamber()-1; // chamber number
	244	Int_t trackIndex = mHit->Track();
	245	if (fNch >= kNCH) {
	246	cerr<<"AliRICHMerger: chamber nr. fNch out of range: "<<fNch<<endl;
	247	cerr<<" track: "<<fTrack<<endl;
	248	continue;
	249	}
	250	iChamber = &(pRICH->Chamber(fNch));
d8a72780	251
b762c2f6	252
	253	//
	254	// If flag is set, create digits only for some particles
	255	//
	256	TParticle current = (TParticle)gAlice->Particle(trackIndex);
	257	if (current->GetPdgCode() >= 50000050)
	258	{
	259	TParticle motherofcurrent = (TParticle)gAlice->Particle(current->GetFirstMother());
	260	particle = motherofcurrent->GetPdgCode();
	261	}
	262	else
	263	{
	264	particle = current->GetPdgCode();
	265	}
d8a72780	266
d8a72780	267
b762c2f6	268	//printf("Flag:%d\n",flag);
	269	//printf("Track:%d\n",track);
	270	//printf("Particle:%d\n",particle);
d8a72780	271
b762c2f6	272	Int_t digitise=1;
d8a72780	273
b762c2f6	274	if (flag == 1)
	275	if(TMath::Abs(particle)==211 \|\| TMath::Abs(particle)==111)
	276	digitise=0;
d8a72780	277
b762c2f6	278	if (flag == 2)
	279	if(TMath::Abs(particle)==321 \|\| TMath::Abs(particle)==130 \|\| TMath::Abs(particle)==310
	280	\|\| TMath::Abs(particle)==311)
	281	digitise=0;
d8a72780	282
b762c2f6	283	if (flag == 3 && TMath::Abs(particle)==2212)
b762c2f6	284	digitise=0;
d8a72780	285
b762c2f6	286	if (flag == 4 && TMath::Abs(particle)==13)
b762c2f6	287	digitise=0;
d8a72780	288
b762c2f6	289	if (flag == 5 && TMath::Abs(particle)==11)
b762c2f6	290	digitise=0;
d8a72780	291
b762c2f6	292	if (flag == 6 && TMath::Abs(particle)==2112)
b762c2f6	293	digitise=0;
d8a72780	294
d8a72780	295
b762c2f6	296	//printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise);
d8a72780	297
d8a72780	298
b762c2f6	299	if (digitise)
b762c2f6	300	{
d8a72780	301
b762c2f6	302	//
	303	// Loop over pad hits
	304	for (AliRICHSDigit* mPad=
	305	(AliRICHSDigit*)pRICH->FirstPad(mHit,pRICH->SDigits());
	306	mPad;
	307	mPad=(AliRICHSDigit*)pRICH->NextPad(pRICH->SDigits()))
	308	{
	309	Int_t iqpad = mPad->QPad(); // charge per pad
	310	fDigits[0]=mPad->PadX();
	311	fDigits[1]=mPad->PadY();
	312	fDigits[2]=iqpad;
	313	fDigits[3]=iqpad;
	314	fDigits[4]=mPad->HitNumber();
d8a72780	315
b762c2f6	316	// build the list of fired pads and update the info
	317	if (Exists(mPad)) {
	318	Update(mPad);
	319	} else {
	320	CreateNew(mPad);
	321	}
	322	} //end loop over clusters
	323	}// track type condition
	324	} // hit loop
	325	} // track loop
d8a72780	326
b762c2f6	327	// open the file with background
d8a72780	328
b762c2f6	329	if (fMerge) {
	330	fSignal = kFALSE;
	331	ntracks = (Int_t)fTrH1->GetEntries();
	332	//
	333	// Loop over tracks
	334	//
	335	for (fTrack = 0; fTrack < ntracks; fTrack++) {
d8a72780	336
b762c2f6	337	if (fHitsBgr) fHitsBgr->Clear();
b762c2f6	338	if (fSDigitsBgr) fSDigitsBgr->Clear();
d8a72780	339
b762c2f6	340	fTrH1->GetEvent(fTrack);
	341	//
	342	// Loop over hits
	343	AliRICHHit* mHit;
	344	for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
	345	{
	346	mHit = (AliRICHHit) (fHitsBgr)[i];
	347	fNch = mHit->Chamber()-1; // chamber number
	348	iChamber = &(pRICH->Chamber(fNch));
	349
d8a72780	350	//
b762c2f6	351	// Loop over pad hits
	352	for (AliRICHSDigit* mPad =
	353	(AliRICHSDigit*)pRICH->FirstPad(mHit,fSDigitsBgr);
	354	mPad;
	355	mPad = (AliRICHSDigit*)pRICH->NextPad(fSDigitsBgr))
	356	{
	357	fDigits[0] = mPad->PadX();
	358	fDigits[1] = mPad->PadY();
	359	fDigits[2] = mPad->QPad();
	360	fDigits[3] = 0;
	361	fDigits[4] = -1; // set hit number to -1 for bgr
d8a72780	362
b762c2f6	363	// build the list of fired pads and update the info
	364	if (!Exists(mPad)) {
	365	CreateNew(mPad);
	366	} else {
	367	Update(mPad);
	368	} // end if !Exists
	369	} //end loop over clusters
	370	} // hit loop
	371	} // track loop
d8a72780	372
b762c2f6	373	TTree *treeK = gAlice->TreeK();
b762c2f6	374	TFile *file = NULL;
d8a72780	375
b762c2f6	376	if (treeK) file = treeK->GetCurrentFile();
	377	file->cd();
	378	} // if fMerge
d8a72780	379
b762c2f6	380	Int_t tracks[kMAXTRACKSPERRICHDIGIT];
	381	Int_t charges[kMAXTRACKSPERRICHDIGIT];
	382	Int_t nentries=fList->GetEntriesFast();
d8a72780	383
b762c2f6	384	for (Int_t nent=0;nent<nentries;nent++) {
	385	AliRICHTransientDigit transDigit=(AliRICHTransientDigit)fList->At(nent);
	386	if (transDigit==0) continue;
	387	Int_t ich=transDigit->GetChamber();
	388	Int_t q=transDigit->Signal();
	389	iChamber=&(pRICH->Chamber(ich));
	390	AliRICHResponse * response=iChamber->GetResponseModel();
	391	Int_t adcmax= (Int_t) response->MaxAdc();
d8a72780	392
d8a72780	393
b762c2f6	394	// add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
	395	//printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
	396	Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());
d8a72780	397
b762c2f6	398	//printf("Pedestal:%d\n",pedestal);
	399	//Int_t pedestal=0;
	400	Float_t treshold = (pedestal + 4*2.2);
d8a72780	401
b762c2f6	402	Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
	403	Float_t noise = gRandom->Gaus(0, meanNoise);
	404	q+=(Int_t)(noise + pedestal);
	405	//q+=(Int_t)(noise);
	406	// magic number to be parametrised !!!
	407	if ( q <= treshold)
	408	{
d8a72780	409	q = q - pedestal;
b762c2f6	410	continue;
	411	}
	412	q = q - pedestal;
	413	if ( q >= adcmax) q=adcmax;
	414	fDigits[0]=transDigit->PadX();
	415	fDigits[1]=transDigit->PadY();
	416	fDigits[2]=q;
	417	fDigits[3]=transDigit->Physics();
	418	fDigits[4]=transDigit->Hit();
	419
	420	Int_t nptracks = transDigit->GetNTracks();
	421
	422	// this was changed to accomodate the real number of tracks
	423	if (nptracks > kMAXTRACKSPERRICHDIGIT) {
	424	printf("Attention - tracks > 10 %d\n",nptracks);
	425	nptracks=kMAXTRACKSPERRICHDIGIT;
	426	}
	427	if (nptracks > 2) {
	428	printf("Attention - tracks > 2 %d \n",nptracks);
	429	//printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
	430	//icat,ich,digits[0],digits[1],q);
	431	}
	432	for (Int_t tr=0;tr<nptracks;tr++) {
	433	tracks[tr]=transDigit->GetTrack(tr);
	434	charges[tr]=transDigit->GetCharge(tr);
	435	} //end loop over list of tracks for one pad
	436	if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
	437	for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
	438	tracks[t]=0;
	439	charges[t]=0;
d8a72780	440	}
b762c2f6	441	}
	442	//write file
	443	//if (ich==2)
	444	//fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
d8a72780	445
b762c2f6	446	// fill digits
	447	pRICH->AddDigits(ich,tracks,charges,fDigits);
	448	}
	449	gAlice->TreeD()->Fill();
	450
	451	fList->Delete();
	452	for(Int_t ii=0;ii<kNCH;++ii) {
	453	if (fHitMap[ii]) {
	454	delete fHitMap[ii];
	455	fHitMap[ii]=0;
d8a72780	456	}
b762c2f6	457	}
d8a72780	458
b762c2f6	459	TClonesArray *richDigits;
	460	for (Int_t k=0;k<kNCH;k++) {
	461	richDigits = pRICH->DigitsAddress(k);
	462	Int_t ndigit=richDigits->GetEntriesFast();
	463	printf ("Chamber %d digits %d \n",k,ndigit);
	464	}
	465	pRICH->ResetDigits(); /// ??? should it be here???
	466	gAlice->TreeD()->Write(0,TObject::kOverwrite);
	467	// reset tree
	468	// gAlice->TreeD()->Reset();
	469	// delete fList; // deleted in dtor
	470	// gObjectTable->Print();
d8a72780	471
d8a72780	472	}