[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 <Riostream.h>
#include <TTree.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TNtuple.h>
#include <TSystem.h>
#include <TROOT.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;
    fHitsBgr     = 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 (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);  this is done in the steering macro
    AliHeader *header = gAlice->GetHeader();
    AliGenEventHeader* mcHeader = header->GenEventHeader();
    fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
    fSpecn  = ((AliGenHijingEventHeader*) mcHeader)->ProjSpectatorsn();
    fSpecp  = ((AliGenHijingEventHeader*) mcHeader)->ProjSpectatorsp();
    // 
    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 = (TFile*)gROOT->GetListOfFiles()->FindObject(fFnBgr);
  if(!file)cerr<<"AliZDCMerger: background file "<<fFnBgr<<" not found\n";
  //    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");

//    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_ROOT/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_ROOT/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] >= 9999) rnd[pl] = 9998;
     //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();
  }
      
}
Commit	Line	Data
	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 <Riostream.h>
	25	#include <TTree.h>
	26	#include <TFile.h>
	27	#include <TDirectory.h>
	28	#include <TNtuple.h>
	29	#include <TSystem.h>
	30	#include <TROOT.h>
	31
	32	// --- AliRoot header files
	33	#include "AliZDCMerger.h"
	34	#include "AliZDC.h"
	35	#include "AliZDCHit.h"
	36	#include "AliZDCMergedHit.h"
	37	#include "AliZDCDigit.h"
	38	#include "AliZDCFragment.h"
	39	#include "AliRun.h"
	40	#include "AliDetector.h"
	41	#include "AliHeader.h"
	42	#include "AliGenEventHeader.h"
	43	#include "AliGenHijingEventHeader.h"
	44
	45	ClassImp(AliZDCMerger)
	46
	47	//int comp(const void i,const void j) {return (int )i - (int )j;}
	48
	49	//____________________________________________________________________________
	50	AliZDCMerger::AliZDCMerger()
	51	{
	52	// Default constructor
	53	//fMerge = kDigitize -> Only digitization
	54	//fMerge = kMerge -> Digitization + Merging
	55	fMerge = kMerge;
	56	fFnBgr = 0;
	57	fBgrFile = 0;
	58	fNEvBgr = 0;
	59	fHitsBgr = 0;
	60	fImpPar = 0;
	61	fSpecn = 0;
	62	fSpecp = 0;
	63	fFreeSpn = 0;
	64	fFreeSpp = 0;
	65	fFnSpecn = 0;
	66	fSpecnFile = 0;
	67	fFnSpecp = 0;
	68	fSpecpFile = 0;
	69	fNMhits = 0;
	70
	71	}
	72
	73	//____________________________________________________________________________
	74	AliZDCMerger::~AliZDCMerger()
	75	{
	76	// Destructor
	77	if (fSpecnFile) delete fSpecnFile;
	78	if (fSpecpFile) delete fSpecpFile;
	79	}
	80
	81	//____________________________________________________________________________
	82	void AliZDCMerger::InitMerging()
	83	{
	84	// Hits tree, impact parameter, num. of spectators n & p
	85	// in background (full Hijing) event
	86	Float_t b;
	87	Background(b, fSpecn, fSpecp);
	88
	89	// Production of nuclear fragments -> num. of FREE spectators n & p
	90	Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
	91
	92	// Extract from spectators distribution the signal events:
	93	// NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
	94	Mixing();
	95	}
	96
	97	//____________________________________________________________________________
	98	void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
	99	{
	100
	101	// --- Open the background file
	102	if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
	103
	104	// --- Read from the TreeE impact parameter (b),
	105	// # of spectators n and p (fSpecn, fSpecp)
	106	fBgrFile->cd();
	107
	108	// // Get AliRun object from file or create it if not on file
	109	// gAlice = (AliRun*)fBgrFile->Get("gAlice");
	110	// if (!gAlice) {
	111	// gAlice = (AliRun*)fBgrFile->Get("gAlice");
	112	// if (gAlice) printf("AliRun object found on file\n");
	113	// if (!gAlice) {
	114	// printf("\n create new gAlice object");
	115	// gAlice = new AliRun("gAlice","Alice test program");
	116	// }
	117	// }
	118
	119	// gAlice->GetEvent(fNEvBgr); this is done in the steering macro
	120	AliHeader *header = gAlice->GetHeader();
	121	AliGenEventHeader* mcHeader = header->GenEventHeader();
	122	fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
	123	fSpecn = ((AliGenHijingEventHeader*) mcHeader)->ProjSpectatorsn();
	124	fSpecp = ((AliGenHijingEventHeader*) mcHeader)->ProjSpectatorsp();
	125	//
	126	printf("\n HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
	127	fNEvBgr,fImpPar,fSpecn,fSpecp);
	128	}
	129
	130	//____________________________________________________________________________
	131	TFile* AliZDCMerger::OpenBgrFile()
	132	{
	133	// Initialise background event
	134	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject(fFnBgr);
	135	if(!file)cerr<<"AliZDCMerger: background file "<<fFnBgr<<" not found\n";
	136	// TFile *file = new TFile(fFnBgr,"UPDATE");
	137	printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
	138	fHitsBgr = new TClonesArray("AliZDCHit",1000);
	139	fMHits = new TClonesArray("AliZDCMergedHit",1000);
	140	return file;
	141	}
	142
	143	//____________________________________________________________________________
	144	void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
	145	Int_t &fFreeSpn, Int_t &fFreeSpp)
	146	{
	147	//printf("\n Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
	148	Int_t j, zz[100], nn[100], nAlpha, ztot, ntot;
	149	AliZDCFragment *frag = new AliZDCFragment(fImpPar);
	150	for(j=0; j<=99; j++){
	151	zz[j] =0;
	152	nn[j] =0;
	153	}
	154
	155	// Fragments generation
	156	frag->GenerateIMF(zz, nAlpha);
	157
	158	// Attach neutrons
	159	ztot=0;
	160	ntot=0;
	161	frag->AttachNeutrons(zz, nn, ztot, ntot);
	162	fFreeSpn = fSpecn-ntot-2*nAlpha;
	163	fFreeSpp = fSpecp-ztot-2*nAlpha;
	164	if(fFreeSpn<0) fFreeSpn=0;
	165	if(fFreeSpp<0) fFreeSpp=0;
	166	//printf("\n 2nAlpha = %d, ztot = %d, ntot = %d\n",2nAlpha, ztot, ntot);
	167	printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
	168	}
	169
	170	//____________________________________________________________________________
	171	void AliZDCMerger::Mixing()
	172	{
	173
	174	//printf("\n AliZDCMerger->Mixing\n");
	175
	176	// ### Background event Hits ###########################################
	177	fBgrFile->cd();
	178	// fBgrFile->ls();
	179
	180	AliZDC zdc = (AliZDC )gAlice->GetModule("ZDC");
	181	// if(zdc) printf("\n Ho trovato lo ZDC!\n");
	182
	183	// fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
	184	// Hits tree
	185	char treeBgrName[20];
	186	sprintf(treeBgrName,"TreeH%d",fNEvBgr);
	187	fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
	188	if(!fTrHBgr){
	189	printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
	190	}
	191	// fTrHBgr->Print();
	192
	193	// Branch address
	194	TBranch *branch;
	195	char branchname[20];
	196	sprintf(branchname,"%s",zdc->GetName());
	197	if(fTrHBgr && fHitsBgr){
	198	// printf("\n fTrHBgr!=0 && fHitsBgr!=0\n");
	199	branch = fTrHBgr->GetBranch(branchname);
	200	if(branch) branch->SetAddress(&fHitsBgr);
	201	}
	202
	203	Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
	204	//printf("\n --- ntracks = %d\n\n", ntracks);
	205
	206	Int_t itrack, nhits, ihit, j, sector[2];
	207	AliZDCHit* zdcHit;
	208	AliZDCMergedHit *mHit;
	209	Float_t mHits[7];
	210	TClonesArray &sdigits = *fMHits; // SDigits TCArray
	211	fNMhits = 0;
	212
	213	// --- Tracks loop
	214	for(itrack=0; itrack<ntracks; itrack++){
	215	fTrHBgr->GetEvent(itrack);
	216
	217	nhits = fHitsBgr->GetEntries();
	218	// printf(" nhits = %d \n", nhits);
	219	for(ihit=0; ihit<nhits; ihit++){
	220	zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
	221
	222	for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
	223	mHits[0] = zdcHit->GetPrimKinEn();
	224	mHits[1] = zdcHit->GetXImpact();
	225	mHits[2] = zdcHit->GetYImpact();
	226	mHits[3] = zdcHit->GetSFlag();
	227	mHits[4] = zdcHit->GetLightPMQ();
	228	mHits[5] = zdcHit->GetLightPMC();
	229	mHits[6] = zdcHit->GetEnergy();
	230	mHit = new AliZDCMergedHit(sector, mHits);
	231	// mHit->Print("");
	232	new((fMHits)[fNMhits]) AliZDCMergedHit(mHit);
	233	new (sdigits[fNMhits]) AliZDCMergedHit(*mHit);
	234	delete mHit;
	235	fNMhits++;
	236	}//Hits loop
	237
	238	} // Tracks loop
	239	//printf(" fNMhits (after bckg) = %d, \n",fNMhits);
	240
	241	// ### Signal event Hits ###########################################
	242	// --- Neutrons
	243	ExtractSignal(1);
	244
	245	// --- Protons
	246	ExtractSignal(2);
	247	//printf(" fNMhits (after signal) = %d \n",fNMhits);
	248
	249	}
	250
	251	//____________________________________________________________________________
	252	void AliZDCMerger::ExtractSignal(Int_t SpecType)
	253	{
	254
	255	// printf("\n Entering in Extract Signal\n");
	256
	257	Int_t numEvents = 0;
	258	if(SpecType == 1){ // --- Signal for spectator neutrons
	259	fFnSpecn = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/ZNsignalntu.root");
	260	fSpecnFile = TFile::Open(fFnSpecn,"R");
	261	fSpecnFile->cd();
	262	printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn);
	263	numEvents = fFreeSpn;
	264	}
	265	else if(SpecType == 2){ // --- Signal for spectator protons
	266	fFnSpecp = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/ZPsignalntu.root");
	267	fSpecpFile = TFile::Open(fFnSpecp,"R");
	268	fSpecpFile->cd();
	269	printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp);
	270	numEvents = fFreeSpp;
	271	}
	272	//printf("\n # of free spectators = %d\n", numEvents);
	273	//printf("\n fNMhits (before adding signal) = %d\n",fNMhits);
	274
	275	TNtuple zdcSignal = (TNtuple) gDirectory->Get("ZDCSignal");
	276	Int_t nentries = (Int_t) zdcSignal->GetEntries();
	277	//printf("\n # entries = %d\n", nentries);
	278
	279	AliZDCMergedHit *mHit;
	280	Float_t *entry, hitsSpec[7];
	281	Int_t pl, i, j, k, iev=0, rnd[125], volume[2];
	282	for(pl=0;pl<125;pl++){
	283	rnd[pl] = 0;
	284	}
	285	for(pl=0;pl<numEvents;pl++){
	286	rnd[pl] = (Int_t) (9999*gRandom->Rndm());
	287	if(rnd[pl] >= 9999) rnd[pl] = 9998;
	288	//printf(" rnd[%d] = %d\n",pl,rnd[pl]);
	289	}
	290	// Sorting vector in ascending order with C function QSORT
	291	qsort((void*)rnd,numEvents,sizeof(Int_t),comp);
	292	//for(pl=0;pl<numEvents;pl++){
	293	////printf(" rnd[%d] = %d\n",pl,rnd[pl]);
	294	//}
	295	do{
	296	for(i=0; i<nentries; i++){
	297	zdcSignal->GetEvent(i);
	298	entry = zdcSignal->GetArgs();
	299	if(entry[0] == rnd[iev]){
	300	for(k=0; k<2; k++) volume[k] = (Int_t) entry[k+1];
	301	for(j=0; j<7; j++){
	302	hitsSpec[j] = entry[j+3];
	303	}
	304	//printf("\n i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
	305	mHit = new AliZDCMergedHit(volume, hitsSpec);
	306	new((fMHits)[fNMhits++]) AliZDCMergedHit(mHit);
	307	delete mHit;
	308	}
	309	else if(entry[0] > rnd[iev]){
	310	iev++;
	311	continue;
	312	}
	313	}
	314	}while(iev<numEvents);
	315
	316	if(SpecType ==1){
	317	//printf("\n fNMhits (after n signal) = %d\n",fNMhits);
	318	fSpecnFile->Close();
	319	}
	320	else if(SpecType == 2){
	321	//printf("\n fNMhits (after p signal) = %d\n",fNMhits);
	322	fSpecpFile->Close();
	323	}
	324
	325	}