///////////////////////////////////////////////////////////////////////////////
// --- ROOT system
-#include <iostream.h>
+#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"
ClassImp(AliZDCMerger)
+//int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
+
//____________________________________________________________________________
AliZDCMerger::AliZDCMerger()
{
// Default constructor
- fMerge = kDigitize;
-// fMerge = kMerge;
+ //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;
fFnSpecp = 0;
fSpecpFile = 0;
fNMhits = 0;
- fMHits = 0;
+
}
//____________________________________________________________________________
AliZDCMerger::~AliZDCMerger()
{
-// Destructor
- if (fBgrFile) delete fBgrFile;
- if (fTrHBgr) delete fTrHBgr;
- if (fTrSDBgr) delete fTrSDBgr;
- if (fHitsBgr) delete fHitsBgr;
+ // Destructor
if (fSpecnFile) delete fSpecnFile;
if (fSpecpFile) delete fSpecpFile;
}
// 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);
+ Background(b, fSpecn, fSpecp);
// Production of nuclear fragments -> num. of FREE spectators n & p
- Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp);
+ Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
// Extract from spectators distribution the signal events:
// NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
{
// --- Open the background file
- if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
+ if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
// --- Read from the TreeE impact parameter (b),
// # of spectators n and p (fSpecn, fSpecp)
// 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)->Spectatorsn();
- fSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
- printf("\n HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
+ 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 = new TFile(fFnBgr,"UPDATE");
+ 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);
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;
+ //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;
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;
+ 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");
+ //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;
+ 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);
// Branch address
TBranch *branch;
char branchname[20];
- sprintf(branchname,"%s",ZDC->GetName());
+ sprintf(branchname,"%s",zdc->GetName());
if(fTrHBgr && fHitsBgr){
// printf("\n fTrHBgr!=0 && fHitsBgr!=0\n");
branch = fTrHBgr->GetBranch(branchname);
}
Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
-// printf("\n --- ntracks = %d\n\n", ntracks);
+ //printf("\n --- ntracks = %d\n\n", ntracks);
Int_t itrack, nhits, ihit, j, sector[2];
AliZDCHit* zdcHit;
- AliZDCMergedHit *MHit;
- Float_t MHits[7];
+ AliZDCMergedHit *mHit;
+ Float_t mHits[7];
TClonesArray &sdigits = *fMHits; // SDigits TCArray
fNMhits = 0;
// --- Tracks loop
for(itrack=0; itrack<ntracks; itrack++){
-// printf(" itrack = %d", itrack);
-// gAlice->ResetHits();
fTrHBgr->GetEvent(itrack);
-// for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
-// zdcHit;
-// zdcHit = (AliZDCHit*)ZDC->NextHit()){
-
nhits = fHitsBgr->GetEntries();
-// nhits = ZDCdp->GetNhits();
// 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;
+ 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("\n fNMhits (after bckg) = %d, \n",fNMhits);
-// fMHits->Dump();
-// AliZDCMergedHit *pippo7 = (AliZDCMergedHit*) fMHits->At(7);
-// pippo7->Dump();
-
+ //printf(" fNMhits (after bckg) = %d, \n",fNMhits);
// ### Signal event Hits ###########################################
// --- Neutrons
// --- Protons
ExtractSignal(2);
-// printf("\n fNMhits (after signal) = %d \n",fNMhits);
+ //printf(" fNMhits (after signal) = %d \n",fNMhits);
}
// printf("\n Entering in Extract Signal\n");
- Int_t NumEvents = 0;
+ 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;
+ 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 --- ExtractSignal x p: file %s opened\n", fFnSpecp);
+ numEvents = fFreeSpp;
}
-// printf("\n # of free spectator = %d\n", NumEvents);
-// printf("\n fNMhits (before adding signal) = %d\n",fNMhits);
+ //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);
- Int_t i, j, k;
+ 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 iev, rnd, Volume[2];
- for(iev=0; iev<NumEvents; iev++){
- rnd = (Int_t) (1000*gRandom->Rndm());
-// printf("\n rnd = %d\n", rnd);
+ 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){
-// printf("\n entry[0] = %f\n", entry[0]);
- for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
+ 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];
+ hitsSpec[j] = entry[j+3];
}
- MHit = new AliZDCMergedHit(Volume, HitsSpec);
- new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
- delete MHit;
+ //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) break;
+ else if(entry[0] > rnd[iev]){
+ iev++;
+ continue;
+ }
}
- }
+ }while(iev<numEvents);
if(SpecType ==1){
-// printf("\n fNMhits (after n signal) = %d\n",fNMhits);
+ //printf("\n fNMhits (after n signal) = %d\n",fNMhits);
fSpecnFile->Close();
}
else if(SpecType == 2){
-// printf("\n fNMhits (after p signal) = %d\n",fNMhits);
+ //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("\n 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);
- 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("\n imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d",
-// 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",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;
- // 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;
-}