[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 <TH1.h>
#include <TSystem.h>

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

ClassImp(AliZDCMerger)

//____________________________________________________________________________
AliZDCMerger::AliZDCMerger()
{
// Default constructor    
    fMerge       = kDigitize;
    fNEvBgr	 = 0;
    fFnBgr       = 0;
    fBgrFile     = 0;
    fTrHBgr      = 0;
    fTrSDBgr     = 0;
    fImpPar      = 0;
    fSpecn       = 0;
    fSpecp       = 0;
    fFreeSpn     = 0;
    fFreeSpp     = 0;
    fFnSpecn     = 0;
    fSpecnFile   = 0;
    fFnSpecp     = 0;
    fSpecpFile   = 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;
    Int_t nspecn, nspecp;
    Background(b, nspecn, nspecp);

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

//____________________________________________________________________________
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 = (AliRun*)fBgrFile->Get("gAlice");
    AliHeader *header = gAlice->GetHeader();
    AliGenEventHeader* mcHeader = header->GenEventHeader();
    fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
    fSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
    fSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
//      printf("\n	HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",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);
    return file;
}

//____________________________________________________________________________
void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
                                 Int_t &fFreeSpn, Int_t &fFreeSpp)
{
    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-Ztot-2*nAlpha;
    fFreeSpp = fSpecp-Ntot-2*nAlpha;
//      printf("\n	Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
}

//____________________________________________________________________________
void AliZDCMerger::Mixing(Int_t fFreeSpn, Int_t fFreeSpp)
{

//      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;    

//    // Read from Event Tree the # of events of the Hijing file 
//    TTree *thead = (TTree *) fBgrFile->Get("TE");
//    if(!te){
//      printf("\n	ERROR! -> File %s does not contain TreeE\n");
//    }
//    fNEvBgr = (Int_t) te->GetEntries();
    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); // TreeH
    if(!fTrHBgr){
      printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
    }    
    Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
//    printf("\n	--- ntracks = %d\n",ntracks);

    // SDigits Tree
    char treeSBgrName[20];
    sprintf(treeSBgrName,"TreeS%d",fNEvBgr);
    fTrSDBgr = (TTree*)gDirectory->Get(treeSBgrName); // TreeH
    if(!fTrSDBgr){
      printf("\n ERROR -> Can't find TreeS%d in background file\n",fNEvBgr);
    }    

    Int_t itrack, i, volume[2], detector, quadrant;
    TClonesArray &sdigits = *fHitsBgr;	// SDigits TCArray
    Float_t hits[10];
    // --- Tracks loop
    for(itrack=0; itrack<ntracks; itrack++){
//       printf("		itrack = %d \n", itrack);
       fTrHBgr->GetEvent(itrack);
       
       fTrack = itrack;
       Int_t NMhits  = 0;
       for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
           zdcHit;
    	   zdcHit = (AliZDCHit*)ZDC->NextHit()){ 

	  for(i=0; i<2; i++) volume[i] = zdcHit->GetVolume(i);
//	  printf("\n		volume[0] = %d volume[1]= %d \n",volume[0], volume[1]);
	  if(volume[0] != 3){
	    // Background hits	
	    hits[7] = zdcHit->GetLightPMQ();
	    hits[8] = zdcHit->GetLightPMC();
//	    printf("\n	Prima ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);

            //Signal hits
	    detector = volume[0];
	    quadrant = volume[1];
	    ExtractSignal(detector, quadrant, fQuadLight, fComLight);	   
	    hits[7] += fQuadLight;
	    hits[8] += fComLight;
//	    printf("\n	Bckg + Signal -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);

            zdcHit->SetLightPMQ(hits[7]);
            zdcHit->SetLightPMC(hits[8]);
	    hits[7] = zdcHit->GetLightPMQ();
	    hits[8] = zdcHit->GetLightPMC();
//	    printf("\n	Dopo ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
	    	    
	    new (sdigits[NMhits++]) AliZDCHit(zdcHit);
//            printf("\n	NMhits = %d\n",NMhits);
	  }
          fBgrFile->cd();
          fTrSDBgr->Fill();
          fTrSDBgr->Write(0,TObject::kOverwrite);
       } //Hits loop
       Digitize();
    } // Tracks loop
        
//    fBgrFile->Close();
        
}

//____________________________________________________________________________
void AliZDCMerger::ExtractSignal(Int_t detector, Int_t quadrant,
                                    Float_t &fQuadLight, Float_t &fComLight)
{

//    printf("\n	Entering ExtractSignal method -> detector = %d quadrant = %d\n",
//           detector, quadrant);

    // Connect spectator n histo's file
    fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignal.root");
    fSpecnFile = TFile::Open(fFnSpecn,"R");
    fSpecnFile->cd();
//    fSpecnFile->ls();
    TH1F *hPMQ1zn = (TH1F*) gDirectory->Get("hPMQ1zn;1");
    TH1F *hPMQ2zn = (TH1F*) gDirectory->Get("hPMQ2zn;1");
    TH1F *hPMQ3zn = (TH1F*) gDirectory->Get("hPMQ3zn;1");
    TH1F *hPMQ4zn = (TH1F*) gDirectory->Get("hPMQ4zn;1");
    TH1F *hPMC1zn = (TH1F*) gDirectory->Get("hPMC1zn;1");
    TH1F *hPMC2zn = (TH1F*) gDirectory->Get("hPMC2zn;1");
    TH1F *hPMC3zn = (TH1F*) gDirectory->Get("hPMC3zn;1");
    TH1F *hPMC4zn = (TH1F*) gDirectory->Get("hPMC4zn;1");
//    Axis_t x = hPMQ1zn -> GetRandom();
//    printf("	hPMQ1zn -> GetRandom() = %f\n",x);
    
    // Connect spectator p histo's file
    fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignal.root");
    fSpecpFile = TFile::Open(fFnSpecp,"R");
    fSpecpFile->cd(); 
//    fSpecpFile->ls();
    TH1F *hPMQ1zp = (TH1F*) gDirectory->Get("hPMQ1zp;1");
    TH1F *hPMQ2zp = (TH1F*) gDirectory->Get("hPMQ2zp;1");
    TH1F *hPMQ3zp = (TH1F*) gDirectory->Get("hPMQ3zp;1");
    TH1F *hPMQ4zp = (TH1F*) gDirectory->Get("hPMQ4zp;1");
    TH1F *hPMC1zp = (TH1F*) gDirectory->Get("hPMC1zp;1");
    TH1F *hPMC2zp = (TH1F*) gDirectory->Get("hPMC2zp;1");
    TH1F *hPMC3zp = (TH1F*) gDirectory->Get("hPMC3zp;1");
    TH1F *hPMC4zp = (TH1F*) gDirectory->Get("hPMC4zp;1");

    if(detector == 1){ // --- ZN
      if(quadrant == 1){
        fQuadLight = (Float_t) hPMQ1zn -> GetRandom();
	fComLight = (Float_t) hPMC1zn -> GetRandom();
      }
      else if(quadrant == 2){
        fQuadLight = (Float_t) hPMQ2zn -> GetRandom();
	fComLight = (Float_t) hPMC2zn -> GetRandom();
      }
      else if(quadrant == 3){
        fQuadLight = (Float_t) hPMQ3zn -> GetRandom();
	fComLight = (Float_t) hPMC3zn -> GetRandom();
      }
      else if(quadrant == 4){
        fQuadLight = (Float_t) hPMQ4zn -> GetRandom();
	fComLight = (Float_t) hPMC4zn -> GetRandom();
      }
    }

    else if(detector == 2){ // --- ZP
      fSpecpFile->cd();	// Connect spectator p histo's file
//      fSpecpFile->ls();
      if(quadrant == 1){
        fQuadLight = (Float_t) hPMQ1zp -> GetRandom();
	fComLight = (Float_t) hPMC1zp -> GetRandom();
      }
      else if(quadrant == 2){
        fQuadLight = (Float_t) hPMQ2zp -> GetRandom();
	fComLight = (Float_t) hPMC2zp -> GetRandom();
      }
      else if(quadrant == 3){
        fQuadLight = (Float_t) hPMQ3zp -> GetRandom();
	fComLight = (Float_t) hPMC3zp -> GetRandom();
      }
      else if(quadrant == 4){
        fQuadLight = (Float_t) hPMQ4zp -> GetRandom();
	fComLight = (Float_t) hPMC4zp -> GetRandom();
      }
    }
//    printf("	---	Exiting ExtractSignal -> fQuadLight = %f, fComLight = %f",
//           fQuadLight,fComLight);

}

//____________________________________________________________________________
void AliZDCMerger::Digitize()
{

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

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

  Int_t itrack, lightQ, lightC, sector[2], digit;
  Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM = 0;
  Int_t i;
  for(i=0; i<4; i++){
     PMQZN[i] = 0;
     PMQZP[i] = 0;
  }

  // ### Digitization "on the flight" (no merging) #######################
  if(fMerge == kDigitize){
//    printf("\n		fMerge == kDigitize\n");
    fTrHBgr = gAlice->TreeH(); // TTree
    TTree *treeH = gAlice->TreeH();
    if(!fTrHBgr){
      printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
    }    
    Int_t ntracks =  (Int_t) treeH->GetEntries();
//    printf("\n	--- ntracks = %d\n",ntracks);
        
    // Loop over tracks
    for(itrack=0; itrack<ntracks; itrack++){
       gAlice->ResetHits();
       gAlice->TreeH()->GetEvent(itrack);
       
       // Loop over hits
       for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
           zdcHit;
    	   zdcHit = (AliZDCHit*)ZDC->NextHit()){ 
	  sector[0] = zdcHit->GetVolume(0);
    	  sector[1] = zdcHit->GetVolume(1);
    	  lightQ    = Int_t(zdcHit->GetLightPMQ());
    	  lightC    = Int_t(zdcHit->GetLightPMC());
//	  printf("\n	Digitise -> DET. = %d, quad = %d", sector[0], sector[1]);
//	  printf("		    PMQ = %d, PMC = %d",  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 
    	    PMZEM = PMZEM + lightC;
    	  }
       } // hits loop
    } // tracks loop
  } // if(fMerge)

  // ### Merging and digitization #######################################
  else if(fMerge == kMerge){
//    printf("\n		fMerge == kMerge\n");

//    Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
//    printf("\n        --- ntracks = %d\n",ntracks);
//        
//    // Loop over tracks
//    for(itrack=0; itrack<ntracks; itrack++){
       fTrHBgr->GetEvent(fTrack);
       
//       printf("\n\n	Track # %d --- Digitise -> ", fTrack);
       // Loop over hits
       for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
           zdcHit;
    	   zdcHit = (AliZDCHit*)ZDC->NextHit()){ 
	  sector[0] = zdcHit->GetVolume(0);
    	  sector[1] = zdcHit->GetVolume(1);
    	  lightQ    = Int_t(zdcHit->GetLightPMQ());
    	  lightC    = Int_t(zdcHit->GetLightPMC());
//	  printf("\n	DET. = %d, quad = %d,PMQ = %d, PMC = %d", 
//	         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 
    	    PMZEM = PMZEM + lightC;
    	  }
       } // hits loop
//    } // tracks loop
  } // if(fMerge)

      
  // 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);
  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	ADCZN[%d] = %d",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	ADCZP = %d",PMCZP,digit);
  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	ADCZP[%d] = %d",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; // Detector = ZZEM
  sector[1] = 0; // Single PM ADC
  digit  = Phe2ADCch(3, 0, PMZEM);
//  printf("\n	ZEM *** PMZEM = %d	ADCZEM = %d",PMZEM,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;
  // Parameters for conversion of light yield in ADC channels
  Float_t fPMGain[3][5];      // PM gain
  Float_t fADCRes;            // ADC conversion factor
  
  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 = 10.;
  
  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>
	28	#include <TH1.h>
53c30603	29	#include <TSystem.h>
359cdddc	30
	31	// --- AliRoot header files
	32	#include "AliZDCMerger.h"
	33	#include "AliZDC.h"
	34	#include "AliZDCHit.h"
	35	#include "AliZDCDigit.h"
	36	#include "AliZDCFragment.h"
	37	#include "AliRun.h"
	38	#include "AliDetector.h"
	39	#include "AliHeader.h"
	40	#include "AliGenEventHeader.h"
53c30603	41	#include "AliGenHijingEventHeader.h"
359cdddc	42
	43	ClassImp(AliZDCMerger)
	44
	45	//____________________________________________________________________________
	46	AliZDCMerger::AliZDCMerger()
	47	{
	48	// Default constructor
	49	fMerge = kDigitize;
	50	fNEvBgr = 0;
	51	fFnBgr = 0;
359cdddc	52	fBgrFile = 0;
	53	fTrHBgr = 0;
	54	fTrSDBgr = 0;
	55	fImpPar = 0;
	56	fSpecn = 0;
	57	fSpecp = 0;
	58	fFreeSpn = 0;
	59	fFreeSpp = 0;
	60	fFnSpecn = 0;
359cdddc	61	fSpecnFile = 0;
359cdddc	62	fFnSpecp = 0;
359cdddc	63	fSpecpFile = 0;
	64
	65	}
	66
	67	//____________________________________________________________________________
	68	AliZDCMerger::~AliZDCMerger()
	69	{
	70	// Destructor
	71	if (fBgrFile) delete fBgrFile;
	72	if (fTrHBgr) delete fTrHBgr;
	73	if (fTrSDBgr) delete fTrSDBgr;
	74	if (fHitsBgr) delete fHitsBgr;
	75	if (fSpecnFile) delete fSpecnFile;
	76	if (fSpecpFile) delete fSpecpFile;
	77	}
	78
	79	//____________________________________________________________________________
	80	void AliZDCMerger::InitMerging()
	81	{
	82	// Hits tree, impact parameter, num. of spectators n & p
	83	// in background (full Hijing) event
	84	Float_t b;
	85	Int_t nspecn, nspecp;
	86	Background(b, nspecn, nspecp);
	87
	88	// Production of nuclear fragments -> num. of FREE spectators n & p
	89	Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp);
	90
	91	// Extract from spectators distribution the signal events:
	92	// NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
	93	Mixing(fFreeSpn, fFreeSpp);
	94	}
	95
	96	//____________________________________________________________________________
	97	void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
	98	{
	99
	100	// --- Open the background file
	101	if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
	102
	103	// --- Read from the TreeE impact parameter (b),
	104	// # of spectators n and p (fSpecn, fSpecp)
	105	fBgrFile->cd();
	106
	107	// Get AliRun object from file or create it if not on file
	108	gAlice = (AliRun*)fBgrFile->Get("gAlice");
	109	if (!gAlice) {
	110	gAlice = (AliRun*)fBgrFile->Get("gAlice");
	111	if (gAlice) printf("AliRun object found on file\n");
	112	if (!gAlice) {
	113	printf("\n create new gAlice object");
	114	gAlice = new AliRun("gAlice","Alice test program");
	115	}
	116	}
	117
	118	// gAlice = (AliRun*)fBgrFile->Get("gAlice");
	119	AliHeader *header = gAlice->GetHeader();
	120	AliGenEventHeader* mcHeader = header->GenEventHeader();
	121	fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
	122	fSpecn = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
	123	fSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
	124	// printf("\n HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
	125	}
	126
127	//____________________________________________________________________________
128	TFile* AliZDCMerger::OpenBgrFile()
129	{
130	// Initialise background event
131	TFile *file = new TFile(fFnBgr,"UPDATE");
132	// printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
133	fHitsBgr = new TClonesArray("AliZDCHit",1000);
134	return file;
135	}
136
137	//____________________________________________________________________________
138	void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
139	Int_t &fFreeSpn, Int_t &fFreeSpp)
140	{
141	Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
142	AliZDCFragment *frag = new AliZDCFragment(fImpPar);
143	for(j=0; j<=99; j++){
144	zz[j] =0;
145	nn[j] =0;
146	}
147
148	// Fragments generation
149	frag->GenerateIMF(zz, nAlpha);
150
151	// Attach neutrons
152	Ztot=0;
153	Ntot=0;
154	frag->AttachNeutrons(zz, nn, Ztot, Ntot);
155	fFreeSpn = fSpecn-Ztot-2*nAlpha;
156	fFreeSpp = fSpecp-Ntot-2*nAlpha;
157	// printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
158	}
159
160	//____________________________________________________________________________
161	void AliZDCMerger::Mixing(Int_t fFreeSpn, Int_t fFreeSpp)
162	{
163
164	// printf("\n AliZDCMerger->Mixing\n");
165
166	// ### Background event Hits ###########################################
167	fBgrFile->cd();
168	// fBgrFile->ls();
169
170	AliZDC ZDC = (AliZDC )gAlice->GetModule("ZDC");
171	// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
172
173	// Hits tree
174	if(fTrHBgr) delete fTrHBgr;
175	fTrHBgr = 0;
176	// SDigits tree
177	if(fTrSDBgr) delete fTrSDBgr;
178	fTrSDBgr = 0;
179
180	// // Read from Event Tree the # of events of the Hijing file
181	// TTree thead = (TTree ) fBgrFile->Get("TE");
182	// if(!te){
183	// printf("\n ERROR! -> File %s does not contain TreeE\n");
184	// }
185	// fNEvBgr = (Int_t) te->GetEntries();
186	fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
187
188	// Hits Tree
189	char treeBgrName[20];
190	sprintf(treeBgrName,"TreeH%d",fNEvBgr);
191	fTrHBgr = (TTree*)gDirectory->Get(treeBgrName); // TreeH
192	if(!fTrHBgr){
193	printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
194	}
195	Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
196	// printf("\n --- ntracks = %d\n",ntracks);
197
198	// SDigits Tree
199	char treeSBgrName[20];
200	sprintf(treeSBgrName,"TreeS%d",fNEvBgr);
201	fTrSDBgr = (TTree*)gDirectory->Get(treeSBgrName); // TreeH
202	if(!fTrSDBgr){
203	printf("\n ERROR -> Can't find TreeS%d in background file\n",fNEvBgr);
204	}
205
206	Int_t itrack, i, volume[2], detector, quadrant;
207	TClonesArray &sdigits = *fHitsBgr; // SDigits TCArray
208	Float_t hits[10];
209	// --- Tracks loop
210	for(itrack=0; itrack<ntracks; itrack++){
211	// printf(" itrack = %d \n", itrack);
212	fTrHBgr->GetEvent(itrack);
213
214	fTrack = itrack;
215	Int_t NMhits = 0;
216	for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
217	zdcHit;
218	zdcHit = (AliZDCHit*)ZDC->NextHit()){
219
220	for(i=0; i<2; i++) volume[i] = zdcHit->GetVolume(i);
221	// printf("\n volume[0] = %d volume[1]= %d \n",volume[0], volume[1]);
222	if(volume[0] != 3){
223	// Background hits
224	hits[7] = zdcHit->GetLightPMQ();
225	hits[8] = zdcHit->GetLightPMC();
226	// printf("\n Prima ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
227
228	//Signal hits
229	detector = volume[0];
230	quadrant = volume[1];
231	ExtractSignal(detector, quadrant, fQuadLight, fComLight);
232	hits[7] += fQuadLight;
233	hits[8] += fComLight;
234	// printf("\n Bckg + Signal -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
235
236	zdcHit->SetLightPMQ(hits[7]);
237	zdcHit->SetLightPMC(hits[8]);
238	hits[7] = zdcHit->GetLightPMQ();
239	hits[8] = zdcHit->GetLightPMC();
240	// printf("\n Dopo ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
241
242	new (sdigits[NMhits++]) AliZDCHit(zdcHit);
243	// printf("\n NMhits = %d\n",NMhits);
244	}
245	fBgrFile->cd();
246	fTrSDBgr->Fill();
247	fTrSDBgr->Write(0,TObject::kOverwrite);
248	} //Hits loop
249	Digitize();
250	} // Tracks loop
251
252	// fBgrFile->Close();
253
254	}
255
256	//____________________________________________________________________________
257	void AliZDCMerger::ExtractSignal(Int_t detector, Int_t quadrant,
258	Float_t &fQuadLight, Float_t &fComLight)
259	{
260
261	// printf("\n Entering ExtractSignal method -> detector = %d quadrant = %d\n",
262	// detector, quadrant);
263
264	// Connect spectator n histo's file
53c30603	265	fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignal.root");
53c30603	266	fSpecnFile = TFile::Open(fFnSpecn,"R");
359cdddc	267	fSpecnFile->cd();
	268	// fSpecnFile->ls();
	269	TH1F hPMQ1zn = (TH1F) gDirectory->Get("hPMQ1zn;1");
	270	TH1F hPMQ2zn = (TH1F) gDirectory->Get("hPMQ2zn;1");
	271	TH1F hPMQ3zn = (TH1F) gDirectory->Get("hPMQ3zn;1");
	272	TH1F hPMQ4zn = (TH1F) gDirectory->Get("hPMQ4zn;1");
	273	TH1F hPMC1zn = (TH1F) gDirectory->Get("hPMC1zn;1");
	274	TH1F hPMC2zn = (TH1F) gDirectory->Get("hPMC2zn;1");
	275	TH1F hPMC3zn = (TH1F) gDirectory->Get("hPMC3zn;1");
	276	TH1F hPMC4zn = (TH1F) gDirectory->Get("hPMC4zn;1");
	277	// Axis_t x = hPMQ1zn -> GetRandom();
	278	// printf(" hPMQ1zn -> GetRandom() = %f\n",x);
	279
	280	// Connect spectator p histo's file
53c30603	281	fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignal.root");
53c30603	282	fSpecpFile = TFile::Open(fFnSpecp,"R");
359cdddc	283	fSpecpFile->cd();
	284	// fSpecpFile->ls();
	285	TH1F hPMQ1zp = (TH1F) gDirectory->Get("hPMQ1zp;1");
	286	TH1F hPMQ2zp = (TH1F) gDirectory->Get("hPMQ2zp;1");
	287	TH1F hPMQ3zp = (TH1F) gDirectory->Get("hPMQ3zp;1");
	288	TH1F hPMQ4zp = (TH1F) gDirectory->Get("hPMQ4zp;1");
	289	TH1F hPMC1zp = (TH1F) gDirectory->Get("hPMC1zp;1");
	290	TH1F hPMC2zp = (TH1F) gDirectory->Get("hPMC2zp;1");
	291	TH1F hPMC3zp = (TH1F) gDirectory->Get("hPMC3zp;1");
	292	TH1F hPMC4zp = (TH1F) gDirectory->Get("hPMC4zp;1");
	293
	294	if(detector == 1){ // --- ZN
	295	if(quadrant == 1){
	296	fQuadLight = (Float_t) hPMQ1zn -> GetRandom();
	297	fComLight = (Float_t) hPMC1zn -> GetRandom();
	298	}
	299	else if(quadrant == 2){
	300	fQuadLight = (Float_t) hPMQ2zn -> GetRandom();
	301	fComLight = (Float_t) hPMC2zn -> GetRandom();
	302	}
	303	else if(quadrant == 3){
	304	fQuadLight = (Float_t) hPMQ3zn -> GetRandom();
	305	fComLight = (Float_t) hPMC3zn -> GetRandom();
	306	}
	307	else if(quadrant == 4){
	308	fQuadLight = (Float_t) hPMQ4zn -> GetRandom();
	309	fComLight = (Float_t) hPMC4zn -> GetRandom();
	310	}
	311	}
	312
	313	else if(detector == 2){ // --- ZP
	314	fSpecpFile->cd(); // Connect spectator p histo's file
	315	// fSpecpFile->ls();
	316	if(quadrant == 1){
	317	fQuadLight = (Float_t) hPMQ1zp -> GetRandom();
	318	fComLight = (Float_t) hPMC1zp -> GetRandom();
	319	}
	320	else if(quadrant == 2){
	321	fQuadLight = (Float_t) hPMQ2zp -> GetRandom();
	322	fComLight = (Float_t) hPMC2zp -> GetRandom();
	323	}
	324	else if(quadrant == 3){
	325	fQuadLight = (Float_t) hPMQ3zp -> GetRandom();
	326	fComLight = (Float_t) hPMC3zp -> GetRandom();
	327	}
	328	else if(quadrant == 4){
	329	fQuadLight = (Float_t) hPMQ4zp -> GetRandom();
	330	fComLight = (Float_t) hPMC4zp -> GetRandom();
	331	}
	332	}
	333	// printf(" --- Exiting ExtractSignal -> fQuadLight = %f, fComLight = %f",
	334	// fQuadLight,fComLight);
	335
	336	}
	337
	338	//____________________________________________________________________________
	339	void AliZDCMerger::Digitize()
	340	{
	341
	342	// printf("\n AliZDCMerger->Digitize()");
	343
	344	AliZDC ZDC = (AliZDC )gAlice->GetModule("ZDC");
	345	// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
	346
347	Int_t itrack, lightQ, lightC, sector[2], digit;
348	Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM = 0;
349	Int_t i;
350	for(i=0; i<4; i++){
351	PMQZN[i] = 0;
352	PMQZP[i] = 0;
353	}
354
355	// ### Digitization "on the flight" (no merging) #######################
356	if(fMerge == kDigitize){
357	// printf("\n fMerge == kDigitize\n");
358	fTrHBgr = gAlice->TreeH(); // TTree
359	TTree *treeH = gAlice->TreeH();
360	if(!fTrHBgr){
361	printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
362	}
363	Int_t ntracks = (Int_t) treeH->GetEntries();
364	// printf("\n --- ntracks = %d\n",ntracks);
365
366	// Loop over tracks
367	for(itrack=0; itrack<ntracks; itrack++){
368	gAlice->ResetHits();
369	gAlice->TreeH()->GetEvent(itrack);
370
371	// Loop over hits
372	for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
373	zdcHit;
374	zdcHit = (AliZDCHit*)ZDC->NextHit()){
375	sector[0] = zdcHit->GetVolume(0);
376	sector[1] = zdcHit->GetVolume(1);
377	lightQ = Int_t(zdcHit->GetLightPMQ());
378	lightC = Int_t(zdcHit->GetLightPMC());
379	// printf("\n Digitise -> DET. = %d, quad = %d", sector[0], sector[1]);
380	// printf(" PMQ = %d, PMC = %d", lightQ, lightC);
381
382	if(sector[0] == 1){ //ZN
383	PMCZN = PMCZN + lightC;
384	PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
385	}
386	else if(sector[0] == 2){ //ZP
387	PMCZP = PMCZP + lightC;
388	PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
389	}
390	else if(sector[0] == 3){ //ZEM
391	PMZEM = PMZEM + lightC;
392	}
393	} // hits loop
394	} // tracks loop
395	} // if(fMerge)
396
397	// ### Merging and digitization #######################################
398	else if(fMerge == kMerge){
399	// printf("\n fMerge == kMerge\n");
400
401	// Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
402	// printf("\n --- ntracks = %d\n",ntracks);
403	//
404	// // Loop over tracks
405	// for(itrack=0; itrack<ntracks; itrack++){
406	fTrHBgr->GetEvent(fTrack);
407
408	// printf("\n\n Track # %d --- Digitise -> ", fTrack);
409	// Loop over hits
410	for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
411	zdcHit;
412	zdcHit = (AliZDCHit*)ZDC->NextHit()){
413	sector[0] = zdcHit->GetVolume(0);
414	sector[1] = zdcHit->GetVolume(1);
415	lightQ = Int_t(zdcHit->GetLightPMQ());
416	lightC = Int_t(zdcHit->GetLightPMC());
417	// printf("\n DET. = %d, quad = %d,PMQ = %d, PMC = %d",
418	// sector[0], sector[1],lightQ, lightC);
419
420	if(sector[0] == 1){ //ZN
421	PMCZN = PMCZN + lightC;
422	PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
423	}
424	else if(sector[0] == 2){ //ZP
425	PMCZP = PMCZP + lightC;
426	PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
427	}
428	else if(sector[0] == 3){ //ZEM
429	PMZEM = PMZEM + lightC;
430	}
431	} // hits loop
432	// } // tracks loop
433	} // if(fMerge)
434
435
436	// Create digits for ZN
437	Int_t PedValue;
438	sector[0] = 1; // Detector = ZN
439	sector[1] = 0; // Common PM ADC
440	digit = Phe2ADCch(1, 0, PMCZN);
441	// printf("\n\n ZN ### PMCZN = %d ADCZN = %d",PMCZN, digit);
442	ZDC->AddDigit(sector, digit);
443	Int_t j;
444	for(j=0; j<4; j++){
445	sector[1] = j+1; // Towers PM ADCs
446	digit = Phe2ADCch(1, j+1, PMQZN[j]);
447	// printf("\n PMQZN[%d] = %d ADCZN[%d] = %d",j,PMQZN[j],j,digit);
448	PedValue = AddPedestal();
449	digit += PedValue;
450	// printf(" PedValue = %d",PedValue);
451	ZDC->AddDigit(sector, digit);
452	}
453	// printf("\n");
454
455	// Create digits for ZP
456	sector[0] = 2; // Detector = ZP
457	sector[1] = 0; // Common PM ADC
458	digit = Phe2ADCch(2, 0, PMCZP);
459	// printf("\n ZP --- PMCZP = %d ADCZP = %d",PMCZP,digit);
460	ZDC->AddDigit(sector, digit);
461	for(j=0; j<4; j++){
462	sector[1] = j+1; // Towers PM ADCs
463	digit = Phe2ADCch(2, j+1, PMQZP[j]);
464	// printf("\n PMQZP[%d] = %d ADCZP[%d] = %d",j,PMQZP[j],j,digit);
465	PedValue = AddPedestal();
466	digit += PedValue;
467	// printf(" PedValue = %d",PedValue);
468	ZDC->AddDigit(sector, digit);
469	}
470	// printf("\n");
471
472	// Create digits for ZEM
473	sector[0] = 3; // Detector = ZZEM
474	sector[1] = 0; // Single PM ADC
475	digit = Phe2ADCch(3, 0, PMZEM);
476	// printf("\n ZEM *** PMZEM = %d ADCZEM = %d",PMZEM,digit);
477	PedValue = AddPedestal();
478	digit += PedValue;
479	// printf(" PedValue = %d\n",PedValue);
480	ZDC->AddDigit(sector, digit);
481
482	}
483
484	//_____________________________________________________________________________
485	Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
486	{
487	// Evaluation of the ADC channel corresponding to the light yield Light
488
489	if(gAlice->GetDebug() > 0){
490	// printf("\n Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
491	}
492
493	Int_t ADCch = 0;
494	// Parameters for conversion of light yield in ADC channels
495	Float_t fPMGain[3][5]; // PM gain
496	Float_t fADCRes; // ADC conversion factor
497
498	Int_t j,i;
499	for(i=0; i<3; i++){
500	for(j=0; j<5; j++){
501	fPMGain[i][j] = 100000.;
502	}
503	}
504	fADCRes = 0.00000064; // ADC Resolution: 250 fC/ADCch
505
506	ADCch = (Int_t) (LightfPMGain[Det-1][Quad]fADCRes);
507
508	return ADCch;
509	}
510
511	//_____________________________________________________________________________
512	Int_t AliZDCMerger::AddPedestal()
513	{
514	// --- Pedestal value -> extracted from a gaussian distribution
515	// obtained from the beam test on the ZEM prototype (Aug. 2000)
516
517	Int_t PedValue;
518	Float_t PedMean = 50.;
519	Float_t PedWidth = 10.;
520
521	PedValue = (Int_t) gRandom->Gaus(PedMean,PedWidth);
522
523	return PedValue;
524	}