[u/mrichter/AliRoot.git] / ZDC / AliZDCMerger.cxx

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


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  			ZDC event merging class                              //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

// --- ROOT system
#include <iostream.h>
#include <TTree.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TNtuple.h>
#include <TSystem.h>

// --- AliRoot header files
#include "AliZDCMerger.h"
#include "AliZDC.h"
#include "AliZDCHit.h"
#include "AliZDCMergedHit.h"
#include "AliZDCDigit.h"
#include "AliZDCFragment.h"
#include "AliRun.h"
#include "AliDetector.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliGenHijingEventHeader.h"

ClassImp(AliZDCMerger)

//int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}

//____________________________________________________________________________
AliZDCMerger::AliZDCMerger()
{
// Default constructor    
    //fMerge       = kDigitize	-> Only digitization
    //fMerge       = kMerge	-> Digitization + Merging
    fMerge       = kMerge;	
    fFnBgr       = 0;
    fBgrFile     = 0;
    fNEvBgr	 = 0;
    fTrHBgr      = 0;
    fHitsBgr     = 0;
    fTrSDBgr     = 0;
    fTrDBgr      = 0;
    fImpPar      = 0;
    fSpecn       = 0;
    fSpecp       = 0;
    fFreeSpn     = 0;
    fFreeSpp     = 0;
    fFnSpecn     = 0;
    fSpecnFile   = 0;
    fFnSpecp     = 0;
    fSpecpFile   = 0;
    fNMhits	 = 0;
     
}

//____________________________________________________________________________
AliZDCMerger::~AliZDCMerger()
{
// Destructor
    if (fBgrFile)    delete fBgrFile;
    if (fTrHBgr)     delete fTrHBgr;
    if (fTrSDBgr)    delete fTrSDBgr;
    if (fHitsBgr)    delete fHitsBgr;
    if (fSpecnFile)  delete fSpecnFile;
    if (fSpecpFile)  delete fSpecpFile;
}

//____________________________________________________________________________
void AliZDCMerger::InitMerging()
{
    // Hits tree, impact parameter, num. of spectators n & p 
    // 		in background (full Hijing) event
    Float_t b;
    Background(b, fSpecn, fSpecp);

    // Production of nuclear fragments -> num. of FREE spectators n & p
    Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
    
    // Extract from spectators distribution the signal events:
    // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p 
    Mixing();
}

//____________________________________________________________________________
void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
{
    
    // --- Open the background file
    if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
    
    // --- Read from the TreeE impact parameter (b),
    //     # of spectators n and p (fSpecn, fSpecp)
    fBgrFile->cd();

//    // Get AliRun object from file or create it if not on file
//    gAlice = (AliRun*)fBgrFile->Get("gAlice");
//    if (!gAlice) {
//      gAlice = (AliRun*)fBgrFile->Get("gAlice");
//      if (gAlice) printf("AliRun object found on file\n");
//      if (!gAlice) {
//	  printf("\n create new gAlice object");
//	  gAlice = new AliRun("gAlice","Alice test program");
//      }
//    }
    
    gAlice->GetEvent(fNEvBgr);
    AliHeader *header = gAlice->GetHeader();
    AliGenEventHeader* mcHeader = header->GenEventHeader();
    fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
    Int_t dSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
    Int_t dSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
    // Until there is only 1 ZDC set the # of spectators must be divided by 2!!!
    fSpecn = dSpecn/2;
    fSpecp = dSpecp/2;
    printf("\n	HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
            fNEvBgr,fImpPar,fSpecn,fSpecp);
}

//____________________________________________________________________________
TFile* AliZDCMerger::OpenBgrFile()
{
    // Initialise background event
    TFile *file = new TFile(fFnBgr,"UPDATE");
    printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
    fHitsBgr = new TClonesArray("AliZDCHit",1000);
    fMHits   = new TClonesArray("AliZDCMergedHit",1000);
    return file;
}

//____________________________________________________________________________
void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
                                 Int_t &fFreeSpn, Int_t &fFreeSpp)
{
    //printf("\n	Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
    Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
    AliZDCFragment *frag = new AliZDCFragment(fImpPar);
    for(j=0; j<=99; j++){
       zz[j] =0;
       nn[j] =0;
    }

    // Fragments generation
    frag->GenerateIMF(zz, nAlpha);

    // Attach neutrons
    Ztot=0;
    Ntot=0;
    frag->AttachNeutrons(zz, nn, Ztot, Ntot);
    fFreeSpn = fSpecn-Ntot-2*nAlpha;
    fFreeSpp = fSpecp-Ztot-2*nAlpha;
    if(fFreeSpn<0) fFreeSpn=0;
    if(fFreeSpp<0) fFreeSpp=0;
    //printf("\n			2*nAlpha = %d, Ztot = %d, Ntot = %d\n",2*nAlpha, Ztot,    Ntot);
    printf("\n	Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
}

//____________________________________________________________________________
void AliZDCMerger::Mixing()
{
    
    //printf("\n	AliZDCMerger->Mixing\n");
       
    // ### Background event Hits ###########################################
    fBgrFile->cd();
//    fBgrFile->ls();
    
   AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
//    if(ZDC) printf("\n	Ho trovato lo ZDC!\n");

    // Hits tree
    if(fTrHBgr) delete fTrHBgr;
    fTrHBgr = 0;        
    // SDigits tree
    if(fTrSDBgr) delete fTrSDBgr;
    fTrSDBgr = 0;    

//    fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
    // Hits tree
    char treeBgrName[20];
    sprintf(treeBgrName,"TreeH%d",fNEvBgr);
    fTrHBgr = (TTree*)gDirectory->Get(treeBgrName); 
    if(!fTrHBgr){
      printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
    }    
//    fTrHBgr->Print();
    
    // Branch address
    TBranch *branch;
    char branchname[20];
    sprintf(branchname,"%s",ZDC->GetName());
    if(fTrHBgr && fHitsBgr){
//      printf("\n	fTrHBgr!=0 && fHitsBgr!=0\n");
      branch = fTrHBgr->GetBranch(branchname);
      if(branch) branch->SetAddress(&fHitsBgr);
    }
    
    Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
    //printf("\n	--- ntracks = %d\n\n", ntracks);
    
    Int_t itrack, nhits, ihit, j, sector[2];
    AliZDCHit* zdcHit;
    AliZDCMergedHit *MHit;
    Float_t MHits[7];
    TClonesArray &sdigits = *fMHits;	// SDigits TCArray
    fNMhits = 0;

    // --- Tracks loop
    for(itrack=0; itrack<ntracks; itrack++){
       fTrHBgr->GetEvent(itrack);
       
       nhits = fHitsBgr->GetEntries();
//       printf(" 	      nhits = %d \n", nhits);
       for(ihit=0; ihit<nhits; ihit++){    
	  zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);

	    for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
	    MHits[0] = zdcHit->GetPrimKinEn();
	    MHits[1] = zdcHit->GetXImpact();
	    MHits[2] = zdcHit->GetYImpact();
	    MHits[3] = zdcHit->GetSFlag();
       	    MHits[4] = zdcHit->GetLightPMQ();
	    MHits[5] = zdcHit->GetLightPMC();
	    MHits[6] = zdcHit->GetEnergy();
	    MHit = new AliZDCMergedHit(sector, MHits);
//	    MHit->Print("");
  	    new((*fMHits)[fNMhits]) AliZDCMergedHit(*MHit);
	    new (sdigits[fNMhits])  AliZDCMergedHit(*MHit);
	    delete MHit;
	    fNMhits++;
	  }//Hits loop
	  
    } // Tracks loop
    //printf("	fNMhits (after bckg) = %d, \n",fNMhits);     
        
    // ### Signal event Hits ###########################################
    // --- Neutrons
    ExtractSignal(1);
        
    // --- Protons
    ExtractSignal(2);
    //printf("	fNMhits (after signal) = %d \n",fNMhits); 
        
}

//____________________________________________________________________________
void AliZDCMerger::ExtractSignal(Int_t SpecType)
{

// printf("\n	Entering in Extract Signal\n");
 
 Int_t NumEvents = 0;
 if(SpecType == 1){		// --- Signal for spectator neutrons
   fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
   fSpecnFile = TFile::Open(fFnSpecn,"R");
   fSpecnFile->cd();
   printf("\n	--- ExtractSignal x n: file %s opened\n", fFnSpecn);
   NumEvents = fFreeSpn;
 }
 else if(SpecType == 2){	// --- Signal for spectator protons
   fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
   fSpecpFile = TFile::Open(fFnSpecp,"R");
   fSpecpFile->cd();
   printf("\n	--- ExtractSignal x p: file %s opened\n", fFnSpecp);
   NumEvents = fFreeSpp;
 }
 //printf("\n		# of free spectators = %d\n", NumEvents);
 //printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);

  TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
  Int_t nentries = (Int_t) ZDCSignal->GetEntries();
  //printf("\n   # entries = %d\n", nentries);
  
  AliZDCMergedHit *MHit; 
  Float_t *entry, HitsSpec[7];
  Int_t pl, i, j, k, iev=0, rnd[125], Volume[2];
  for(pl=0;pl<125;pl++){
     rnd[pl] = 0;
  }
  for(pl=0;pl<NumEvents;pl++){
     rnd[pl] = (Int_t) (9999*gRandom->Rndm());
     if(rnd[pl] >= 9998) rnd[pl] = 9997;
     //printf("	rnd[%d] = %d\n",pl,rnd[pl]);     
  }
  // Sorting vector in ascending order with C function QSORT 
  qsort((void*)rnd,NumEvents,sizeof(Int_t),comp);
  //for(pl=0;pl<NumEvents;pl++){
  ////printf("	rnd[%d] = %d\n",pl,rnd[pl]);     
  //}
  do{
     for(i=0; i<nentries; i++){  
  	ZDCSignal->GetEvent(i);
  	entry = ZDCSignal->GetArgs();
  	if(entry[0] == rnd[iev]){
          for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
          for(j=0; j<7; j++){
             HitsSpec[j] = entry[j+3];
          }
          //printf("\n   i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
	  MHit = new AliZDCMergedHit(Volume, HitsSpec);
	  new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
	  delete MHit;
  	}
  	else if(entry[0] > rnd[iev]){
	  iev++;
	  continue;
	}
     }
  }while(iev<NumEvents);
  
  if(SpecType ==1){
    //printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
    fSpecnFile->Close();
  }
  else if(SpecType == 2){
    //printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
    fSpecpFile->Close();
  }
      
}

//____________________________________________________________________________
void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
{

  printf("\n	AliZDCMerger->Digitize()");

  AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
//  if(ZDC) printf("\n 	Ho trovato lo ZDC!\n");

  Int_t lightQ, lightC, sector[2], digit;
  Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
  Int_t i;
  for(i=0; i<4; i++){
     PMQZN[i] = 0;
     PMQZP[i] = 0;
  }
  
  AliZDCMergedHit *MHit;
  Int_t imhit;
  printf("	   fNMHits = %d\n", fNMhits);
  // Loop over SDigits
  for(imhit=0; imhit<fNMhits; imhit++){
     
     MHit = (AliZDCMergedHit*) fMHits->UncheckedAt(imhit);
     sector[0] = MHit->GetSector(0);
     sector[1] = MHit->GetSector(1);
   // tmp -> c'erano i quadranti cannati!
  if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
     printf("\n  *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
 	      sector[0], sector[1]);
       sector[1] = 0;
   }
     lightQ    = Int_t(MHit->GetLightPMQ());
     lightC    = Int_t(MHit->GetLightPMC());
//     printf("    imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n", 
//      imhit,sector[0], sector[1],lightQ, lightC);
     
     if(sector[0] == 1){   	   //ZN 
       PMCZN = PMCZN + lightC;
       PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
     }
     else if(sector[0] == 2){	   //ZP 
       PMCZP = PMCZP + lightC;
       PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
     }
     else if(sector[0] == 3){	   //ZEM 
       if(sector[1] ==1) PMZEM1 = PMZEM1 + lightC;
       else              PMZEM2 = PMZEM2 + lightQ;
     }
  } // SDigits loop
      
  // ### Digits creation ###############################################
  // Create digits for ZN
  Int_t PedValue;
  sector[0] = 1; // Detector = ZN
  sector[1] = 0; // Common PM ADC
  digit = Phe2ADCch(1, 0, PMCZN);
  printf("\n\n	ZN ###	PMCZN = %d	ADCZN = %d",PMCZN, digit);
  PedValue = AddPedestal();
  digit += PedValue;
//  printf("	PedValue = %d",PedValue);
  ZDC->AddDigit(sector, digit);
  Int_t j;
  for(j=0; j<4; j++){
    sector[1] = j+1; // Towers PM ADCs
    digit = Phe2ADCch(1, j+1, PMQZN[j]);
    printf("\n		PMQZN[%d] = %d	phe	ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
    PedValue = AddPedestal();
    digit += PedValue;
//    printf("	PedValue = %d",PedValue);
    ZDC->AddDigit(sector, digit);
  }
  printf("\n");
  
  // Create digits for ZP
  sector[0] = 2; // Detector = ZP
  sector[1] = 0; // Common PM ADC
  digit = Phe2ADCch(2, 0, PMCZP);
  printf("\n	ZP --- PMCZP = %d	phe	ADCZP = %d ADCch",PMCZP,digit);
  PedValue = AddPedestal();
  digit += PedValue;
//  printf("	PedValue = %d",PedValue);
  ZDC->AddDigit(sector, digit);
  for(j=0; j<4; j++){
    sector[1] = j+1; // Towers PM ADCs
    digit = Phe2ADCch(2, j+1, PMQZP[j]);
    printf("\n	       PMQZP[%d] = %d	phe	ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
    PedValue = AddPedestal();
    digit += PedValue;
//    printf("	PedValue = %d",PedValue);
    ZDC->AddDigit(sector, digit);
  }
  printf("\n");
  
  // Create digits for ZEM
  sector[0] = 3; 
  sector[1] = 1; // Detector = ZEM1
  digit  = Phe2ADCch(3, 1, PMZEM1);
  printf("\n	ZEM *** PMZEM1 = %d      phe     ADCZEM1 = %d ADCch",PMZEM1,digit);
  PedValue = AddPedestal();
  digit += PedValue;
//  printf("  PedValue = %d\n",PedValue);
  ZDC->AddDigit(sector, digit); 
  sector[1] = 2; // Detector = ZEM2
  digit  = Phe2ADCch(3, 2, PMZEM2);
  printf("\n	ZEM *** PMZEM2 = %d      phe     ADCZEM2 = %d ADCch\n",PMZEM2,digit);
  PedValue = AddPedestal();
  digit += PedValue;
//  printf("  PedValue = %d\n",PedValue);
  ZDC->AddDigit(sector, digit); 
}

//_____________________________________________________________________________
Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
{
  // Evaluation of the ADC channel corresponding to the light yield Light

  if(gAlice->GetDebug() > 0){
//    printf("\n  Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
  }
  
  Int_t ADCch = 0;
  
  Int_t j,i;
  for(i=0; i<3; i++){
     for(j=0; j<5; j++){
        fPMGain[i][j]   = 100000.;
     }
  }
  fADCRes   = 0.00000064; // ADC Resolution: 250 fC/ADCch
  
  ADCch = (Int_t) (Light*fPMGain[Det-1][Quad]*fADCRes);
     
  return ADCch;
}

//_____________________________________________________________________________
Int_t AliZDCMerger::AddPedestal()
{
  // --- Pedestal value -> extracted from a gaussian distribution
  // obtained from the beam test on the ZEM prototype (Aug. 2000)
  
  Int_t PedValue;
  Float_t PedMean  = 50.;
  Float_t PedWidth = 5.;
  
  PedValue    = (Int_t) gRandom->Gaus(PedMean,PedWidth);
  
  return PedValue;
}
Commit	Line	Data
359cdddc	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// ZDC event merging class //
	20	// //
	21	///////////////////////////////////////////////////////////////////////////////
	22
	23	// --- ROOT system
	24	#include <iostream.h>
	25	#include <TTree.h>
	26	#include <TFile.h>
	27	#include <TDirectory.h>
1acbaaa8	28	#include <TNtuple.h>
53c30603	29	#include <TSystem.h>
359cdddc	30
	31	// --- AliRoot header files
	32	#include "AliZDCMerger.h"
	33	#include "AliZDC.h"
	34	#include "AliZDCHit.h"
1acbaaa8	35	#include "AliZDCMergedHit.h"
359cdddc	36	#include "AliZDCDigit.h"
	37	#include "AliZDCFragment.h"
	38	#include "AliRun.h"
	39	#include "AliDetector.h"
	40	#include "AliHeader.h"
	41	#include "AliGenEventHeader.h"
53c30603	42	#include "AliGenHijingEventHeader.h"
359cdddc	43
	44	ClassImp(AliZDCMerger)
	45
1450a7cd	46	//int comp(const void i,const void j) {return (int )i - (int )j;}
1450a7cd	47
359cdddc	48	//____________________________________________________________________________
	49	AliZDCMerger::AliZDCMerger()
	50	{
	51	// Default constructor
1450a7cd	52	//fMerge = kDigitize -> Only digitization
	53	//fMerge = kMerge -> Digitization + Merging
	54	fMerge = kMerge;
359cdddc	55	fFnBgr = 0;
359cdddc	56	fBgrFile = 0;
1acbaaa8	57	fNEvBgr = 0;
359cdddc	58	fTrHBgr = 0;
2685bf00	59	fHitsBgr = 0;
359cdddc	60	fTrSDBgr = 0;
2685bf00	61	fTrDBgr = 0;
359cdddc	62	fImpPar = 0;
	63	fSpecn = 0;
	64	fSpecp = 0;
	65	fFreeSpn = 0;
	66	fFreeSpp = 0;
	67	fFnSpecn = 0;
359cdddc	68	fSpecnFile = 0;
359cdddc	69	fFnSpecp = 0;
359cdddc	70	fSpecpFile = 0;
1acbaaa8	71	fNMhits = 0;
1450a7cd	72
359cdddc	73	}
	74
	75	//____________________________________________________________________________
	76	AliZDCMerger::~AliZDCMerger()
	77	{
	78	// Destructor
	79	if (fBgrFile) delete fBgrFile;
	80	if (fTrHBgr) delete fTrHBgr;
	81	if (fTrSDBgr) delete fTrSDBgr;
	82	if (fHitsBgr) delete fHitsBgr;
	83	if (fSpecnFile) delete fSpecnFile;
	84	if (fSpecpFile) delete fSpecpFile;
	85	}
	86
	87	//____________________________________________________________________________
	88	void AliZDCMerger::InitMerging()
	89	{
	90	// Hits tree, impact parameter, num. of spectators n & p
	91	// in background (full Hijing) event
	92	Float_t b;
1450a7cd	93	Background(b, fSpecn, fSpecp);
359cdddc	94
359cdddc	95	// Production of nuclear fragments -> num. of FREE spectators n & p
1450a7cd	96	Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
359cdddc	97
	98	// Extract from spectators distribution the signal events:
	99	// NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
1acbaaa8	100	Mixing();
359cdddc	101	}
	102
	103	//____________________________________________________________________________
	104	void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
	105	{
	106
	107	// --- Open the background file
	108	if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
	109
	110	// --- Read from the TreeE impact parameter (b),
	111	// # of spectators n and p (fSpecn, fSpecp)
	112	fBgrFile->cd();
	113
1acbaaa8	114	// // Get AliRun object from file or create it if not on file
359cdddc	115	// gAlice = (AliRun*)fBgrFile->Get("gAlice");
1acbaaa8	116	// if (!gAlice) {
	117	// gAlice = (AliRun*)fBgrFile->Get("gAlice");
	118	// if (gAlice) printf("AliRun object found on file\n");
	119	// if (!gAlice) {
	120	// printf("\n create new gAlice object");
	121	// gAlice = new AliRun("gAlice","Alice test program");
	122	// }
	123	// }
1450a7cd	124
1450a7cd	125	gAlice->GetEvent(fNEvBgr);
359cdddc	126	AliHeader *header = gAlice->GetHeader();
	127	AliGenEventHeader* mcHeader = header->GenEventHeader();
	128	fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
1450a7cd	129	Int_t dSpecn = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
	130	Int_t dSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
	131	// Until there is only 1 ZDC set the # of spectators must be divided by 2!!!
	132	fSpecn = dSpecn/2;
	133	fSpecp = dSpecp/2;
	134	printf("\n HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
	135	fNEvBgr,fImpPar,fSpecn,fSpecp);
359cdddc	136	}
	137
	138	//____________________________________________________________________________
	139	TFile* AliZDCMerger::OpenBgrFile()
	140	{
	141	// Initialise background event
	142	TFile *file = new TFile(fFnBgr,"UPDATE");
1acbaaa8	143	printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
359cdddc	144	fHitsBgr = new TClonesArray("AliZDCHit",1000);
1acbaaa8	145	fMHits = new TClonesArray("AliZDCMergedHit",1000);
359cdddc	146	return file;
	147	}
	148
	149	//____________________________________________________________________________
	150	void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
	151	Int_t &fFreeSpn, Int_t &fFreeSpp)
	152	{
1450a7cd	153	//printf("\n Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
359cdddc	154	Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
	155	AliZDCFragment *frag = new AliZDCFragment(fImpPar);
	156	for(j=0; j<=99; j++){
	157	zz[j] =0;
	158	nn[j] =0;
	159	}
	160
	161	// Fragments generation
	162	frag->GenerateIMF(zz, nAlpha);
	163
	164	// Attach neutrons
	165	Ztot=0;
	166	Ntot=0;
	167	frag->AttachNeutrons(zz, nn, Ztot, Ntot);
1450a7cd	168	fFreeSpn = fSpecn-Ntot-2*nAlpha;
	169	fFreeSpp = fSpecp-Ztot-2*nAlpha;
	170	if(fFreeSpn<0) fFreeSpn=0;
	171	if(fFreeSpp<0) fFreeSpp=0;
	172	//printf("\n 2nAlpha = %d, Ztot = %d, Ntot = %d\n",2nAlpha, Ztot, Ntot);
1acbaaa8	173	printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
359cdddc	174	}
	175
	176	//____________________________________________________________________________
1acbaaa8	177	void AliZDCMerger::Mixing()
359cdddc	178	{
1acbaaa8	179
1450a7cd	180	//printf("\n AliZDCMerger->Mixing\n");
1acbaaa8	181
359cdddc	182	// ### Background event Hits ###########################################
	183	fBgrFile->cd();
	184	// fBgrFile->ls();
	185
1acbaaa8	186	AliZDC ZDC = (AliZDC )gAlice->GetModule("ZDC");
359cdddc	187	// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
	188
	189	// Hits tree
	190	if(fTrHBgr) delete fTrHBgr;
1acbaaa8	191	fTrHBgr = 0;
359cdddc	192	// SDigits tree
	193	if(fTrSDBgr) delete fTrSDBgr;
	194	fTrSDBgr = 0;
	195
1acbaaa8	196	// fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
1acbaaa8	197	// Hits tree
359cdddc	198	char treeBgrName[20];
359cdddc	199	sprintf(treeBgrName,"TreeH%d",fNEvBgr);
1acbaaa8	200	fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
359cdddc	201	if(!fTrHBgr){
	202	printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
	203	}
1acbaaa8	204	// fTrHBgr->Print();
	205
	206	// Branch address
	207	TBranch *branch;
	208	char branchname[20];
	209	sprintf(branchname,"%s",ZDC->GetName());
	210	if(fTrHBgr && fHitsBgr){
	211	// printf("\n fTrHBgr!=0 && fHitsBgr!=0\n");
	212	branch = fTrHBgr->GetBranch(branchname);
	213	if(branch) branch->SetAddress(&fHitsBgr);
	214	}
	215
359cdddc	216	Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
1450a7cd	217	//printf("\n --- ntracks = %d\n\n", ntracks);
1acbaaa8	218
	219	Int_t itrack, nhits, ihit, j, sector[2];
	220	AliZDCHit* zdcHit;
	221	AliZDCMergedHit *MHit;
	222	Float_t MHits[7];
	223	TClonesArray &sdigits = *fMHits; // SDigits TCArray
	224	fNMhits = 0;
359cdddc	225
359cdddc	226	// --- Tracks loop
359cdddc	227	for(itrack=0; itrack<ntracks; itrack++){
359cdddc	228	fTrHBgr->GetEvent(itrack);
359cdddc	229
1acbaaa8	230	nhits = fHitsBgr->GetEntries();
1acbaaa8	231	// printf(" nhits = %d \n", nhits);
	232	for(ihit=0; ihit<nhits; ihit++){
	233	zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
	234
	235	for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
	236	MHits[0] = zdcHit->GetPrimKinEn();
	237	MHits[1] = zdcHit->GetXImpact();
	238	MHits[2] = zdcHit->GetYImpact();
	239	MHits[3] = zdcHit->GetSFlag();
1450a7cd	240	MHits[4] = zdcHit->GetLightPMQ();
1acbaaa8	241	MHits[5] = zdcHit->GetLightPMC();
	242	MHits[6] = zdcHit->GetEnergy();
	243	MHit = new AliZDCMergedHit(sector, MHits);
	244	// MHit->Print("");
	245	new((fMHits)[fNMhits]) AliZDCMergedHit(MHit);
	246	new (sdigits[fNMhits]) AliZDCMergedHit(*MHit);
	247	delete MHit;
	248	fNMhits++;
	249	}//Hits loop
	250
359cdddc	251	} // Tracks loop
1450a7cd	252	//printf(" fNMhits (after bckg) = %d, \n",fNMhits);
1acbaaa8	253
	254	// ### Signal event Hits ###########################################
	255	// --- Neutrons
	256	ExtractSignal(1);
359cdddc	257
1acbaaa8	258	// --- Protons
1acbaaa8	259	ExtractSignal(2);
1450a7cd	260	//printf(" fNMhits (after signal) = %d \n",fNMhits);
359cdddc	261
	262	}
	263
	264	//____________________________________________________________________________
1acbaaa8	265	void AliZDCMerger::ExtractSignal(Int_t SpecType)
359cdddc	266	{
359cdddc	267
1acbaaa8	268	// printf("\n Entering in Extract Signal\n");
	269
	270	Int_t NumEvents = 0;
	271	if(SpecType == 1){ // --- Signal for spectator neutrons
	272	fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
	273	fSpecnFile = TFile::Open(fFnSpecn,"R");
	274	fSpecnFile->cd();
1450a7cd	275	printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn);
1acbaaa8	276	NumEvents = fFreeSpn;
	277	}
	278	else if(SpecType == 2){ // --- Signal for spectator protons
	279	fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
	280	fSpecpFile = TFile::Open(fFnSpecp,"R");
	281	fSpecpFile->cd();
1450a7cd	282	printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp);
1acbaaa8	283	NumEvents = fFreeSpp;
1acbaaa8	284	}
1450a7cd	285	//printf("\n # of free spectators = %d\n", NumEvents);
1450a7cd	286	//printf("\n fNMhits (before adding signal) = %d\n",fNMhits);
1acbaaa8	287
	288	TNtuple ZDCSignal = (TNtuple) gDirectory->Get("ZDCSignal");
	289	Int_t nentries = (Int_t) ZDCSignal->GetEntries();
1450a7cd	290	//printf("\n # entries = %d\n", nentries);
1acbaaa8	291
	292	AliZDCMergedHit *MHit;
	293	Float_t *entry, HitsSpec[7];
1450a7cd	294	Int_t pl, i, j, k, iev=0, rnd[125], Volume[2];
	295	for(pl=0;pl<125;pl++){
	296	rnd[pl] = 0;
	297	}
	298	for(pl=0;pl<NumEvents;pl++){
	299	rnd[pl] = (Int_t) (9999*gRandom->Rndm());
	300	if(rnd[pl] >= 9998) rnd[pl] = 9997;
	301	//printf(" rnd[%d] = %d\n",pl,rnd[pl]);
	302	}
	303	// Sorting vector in ascending order with C function QSORT
	304	qsort((void*)rnd,NumEvents,sizeof(Int_t),comp);
	305	//for(pl=0;pl<NumEvents;pl++){
	306	////printf(" rnd[%d] = %d\n",pl,rnd[pl]);
	307	//}
	308	do{
1acbaaa8	309	for(i=0; i<nentries; i++){
	310	ZDCSignal->GetEvent(i);
	311	entry = ZDCSignal->GetArgs();
1450a7cd	312	if(entry[0] == rnd[iev]){
1acbaaa8	313	for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
	314	for(j=0; j<7; j++){
	315	HitsSpec[j] = entry[j+3];
	316	}
1450a7cd	317	//printf("\n i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
1acbaaa8	318	MHit = new AliZDCMergedHit(Volume, HitsSpec);
	319	new((fMHits)[fNMhits++]) AliZDCMergedHit(MHit);
	320	delete MHit;
	321	}
1450a7cd	322	else if(entry[0] > rnd[iev]){
	323	iev++;
	324	continue;
	325	}
1acbaaa8	326	}
1450a7cd	327	}while(iev<NumEvents);
1acbaaa8	328
1acbaaa8	329	if(SpecType ==1){
1450a7cd	330	//printf("\n fNMhits (after n signal) = %d\n",fNMhits);
1acbaaa8	331	fSpecnFile->Close();
	332	}
	333	else if(SpecType == 2){
1450a7cd	334	//printf("\n fNMhits (after p signal) = %d\n",fNMhits);
1acbaaa8	335	fSpecpFile->Close();
	336	}
	337
359cdddc	338	}
	339
	340	//____________________________________________________________________________
1acbaaa8	341	void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
359cdddc	342	{
359cdddc	343
1450a7cd	344	printf("\n AliZDCMerger->Digitize()");
359cdddc	345
	346	AliZDC ZDC = (AliZDC )gAlice->GetModule("ZDC");
	347	// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
	348
1acbaaa8	349	Int_t lightQ, lightC, sector[2], digit;
1acbaaa8	350	Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
359cdddc	351	Int_t i;
	352	for(i=0; i<4; i++){
	353	PMQZN[i] = 0;
	354	PMQZP[i] = 0;
	355	}
1acbaaa8	356
	357	AliZDCMergedHit *MHit;
	358	Int_t imhit;
1450a7cd	359	printf(" fNMHits = %d\n", fNMhits);
1acbaaa8	360	// Loop over SDigits
	361	for(imhit=0; imhit<fNMhits; imhit++){
	362
	363	MHit = (AliZDCMergedHit*) fMHits->UncheckedAt(imhit);
	364	sector[0] = MHit->GetSector(0);
	365	sector[1] = MHit->GetSector(1);
1450a7cd	366	// tmp -> c'erano i quadranti cannati!
1acbaaa8	367	if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
1450a7cd	368	printf("\n *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
1450a7cd	369	sector[0], sector[1]);
1acbaaa8	370	sector[1] = 0;
	371	}
	372	lightQ = Int_t(MHit->GetLightPMQ());
	373	lightC = Int_t(MHit->GetLightPMC());
1450a7cd	374	// printf(" imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n",
1acbaaa8	375	// imhit,sector[0], sector[1],lightQ, lightC);
	376
	377	if(sector[0] == 1){ //ZN
	378	PMCZN = PMCZN + lightC;
	379	PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
	380	}
	381	else if(sector[0] == 2){ //ZP
	382	PMCZP = PMCZP + lightC;
	383	PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
	384	}
	385	else if(sector[0] == 3){ //ZEM
	386	if(sector[1] ==1) PMZEM1 = PMZEM1 + lightC;
	387	else PMZEM2 = PMZEM2 + lightQ;
	388	}
	389	} // SDigits loop
359cdddc	390
1acbaaa8	391	// ### Digits creation ###############################################
359cdddc	392	// Create digits for ZN
	393	Int_t PedValue;
	394	sector[0] = 1; // Detector = ZN
	395	sector[1] = 0; // Common PM ADC
	396	digit = Phe2ADCch(1, 0, PMCZN);
1450a7cd	397	printf("\n\n ZN ### PMCZN = %d ADCZN = %d",PMCZN, digit);
1acbaaa8	398	PedValue = AddPedestal();
	399	digit += PedValue;
	400	// printf(" PedValue = %d",PedValue);
359cdddc	401	ZDC->AddDigit(sector, digit);
	402	Int_t j;
	403	for(j=0; j<4; j++){
	404	sector[1] = j+1; // Towers PM ADCs
	405	digit = Phe2ADCch(1, j+1, PMQZN[j]);
1450a7cd	406	printf("\n PMQZN[%d] = %d phe ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
359cdddc	407	PedValue = AddPedestal();
	408	digit += PedValue;
	409	// printf(" PedValue = %d",PedValue);
	410	ZDC->AddDigit(sector, digit);
	411	}
1450a7cd	412	printf("\n");
359cdddc	413
	414	// Create digits for ZP
	415	sector[0] = 2; // Detector = ZP
	416	sector[1] = 0; // Common PM ADC
	417	digit = Phe2ADCch(2, 0, PMCZP);
1450a7cd	418	printf("\n ZP --- PMCZP = %d phe ADCZP = %d ADCch",PMCZP,digit);
1acbaaa8	419	PedValue = AddPedestal();
1acbaaa8	420	digit += PedValue;
1450a7cd	421	// printf(" PedValue = %d",PedValue);
359cdddc	422	ZDC->AddDigit(sector, digit);
	423	for(j=0; j<4; j++){
	424	sector[1] = j+1; // Towers PM ADCs
	425	digit = Phe2ADCch(2, j+1, PMQZP[j]);
1450a7cd	426	printf("\n PMQZP[%d] = %d phe ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
359cdddc	427	PedValue = AddPedestal();
	428	digit += PedValue;
	429	// printf(" PedValue = %d",PedValue);
	430	ZDC->AddDigit(sector, digit);
	431	}
1450a7cd	432	printf("\n");
359cdddc	433
359cdddc	434	// Create digits for ZEM
1acbaaa8	435	sector[0] = 3;
	436	sector[1] = 1; // Detector = ZEM1
	437	digit = Phe2ADCch(3, 1, PMZEM1);
1450a7cd	438	printf("\n ZEM *** PMZEM1 = %d phe ADCZEM1 = %d ADCch",PMZEM1,digit);
1acbaaa8	439	PedValue = AddPedestal();
	440	digit += PedValue;
	441	// printf(" PedValue = %d\n",PedValue);
	442	ZDC->AddDigit(sector, digit);
	443	sector[1] = 2; // Detector = ZEM2
	444	digit = Phe2ADCch(3, 2, PMZEM2);
1450a7cd	445	printf("\n ZEM *** PMZEM2 = %d phe ADCZEM2 = %d ADCch\n",PMZEM2,digit);
1acbaaa8	446	PedValue = AddPedestal();
	447	digit += PedValue;
	448	// printf(" PedValue = %d\n",PedValue);
359cdddc	449	ZDC->AddDigit(sector, digit);
359cdddc	450	}
	451
	452	//_____________________________________________________________________________
	453	Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
	454	{
	455	// Evaluation of the ADC channel corresponding to the light yield Light
	456
	457	if(gAlice->GetDebug() > 0){
	458	// printf("\n Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
	459	}
	460
	461	Int_t ADCch = 0;
359cdddc	462
	463	Int_t j,i;
	464	for(i=0; i<3; i++){
	465	for(j=0; j<5; j++){
	466	fPMGain[i][j] = 100000.;
	467	}
	468	}
	469	fADCRes = 0.00000064; // ADC Resolution: 250 fC/ADCch
	470
	471	ADCch = (Int_t) (LightfPMGain[Det-1][Quad]fADCRes);
	472
	473	return ADCch;
	474	}
	475
	476	//_____________________________________________________________________________
	477	Int_t AliZDCMerger::AddPedestal()
	478	{
	479	// --- Pedestal value -> extracted from a gaussian distribution
	480	// obtained from the beam test on the ZEM prototype (Aug. 2000)
	481
	482	Int_t PedValue;
	483	Float_t PedMean = 50.;
1450a7cd	484	Float_t PedWidth = 5.;
359cdddc	485
	486	PedValue = (Int_t) gRandom->Gaus(PedMean,PedWidth);
	487
	488	return PedValue;
	489	}