1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // ZDC event merging class //
21 ///////////////////////////////////////////////////////////////////////////////
24 #include <Riostream.h>
27 #include <TDirectory.h>
32 // --- AliRoot header files
33 #include "AliZDCMerger.h"
35 #include "AliZDCHit.h"
36 #include "AliZDCMergedHit.h"
37 #include "AliZDCDigit.h"
38 #include "AliZDCFragment.h"
40 #include "AliDetector.h"
41 #include "AliHeader.h"
42 #include "AliGenEventHeader.h"
43 #include "AliGenHijingEventHeader.h"
45 ClassImp(AliZDCMerger)
47 //int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
49 //____________________________________________________________________________
50 AliZDCMerger::AliZDCMerger()
52 // Default constructor
53 //fMerge = kDigitize -> Only digitization
54 //fMerge = kMerge -> Digitization + Merging
73 //____________________________________________________________________________
74 AliZDCMerger::~AliZDCMerger()
77 if (fSpecnFile) delete fSpecnFile;
78 if (fSpecpFile) delete fSpecpFile;
81 //____________________________________________________________________________
82 void AliZDCMerger::InitMerging()
84 // Hits tree, impact parameter, num. of spectators n & p
85 // in background (full Hijing) event
87 Background(b, fSpecn, fSpecp);
89 // Production of nuclear fragments -> num. of FREE spectators n & p
90 Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
92 // Extract from spectators distribution the signal events:
93 // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
97 //____________________________________________________________________________
98 void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
101 // --- Open the background file
102 if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
104 // --- Read from the TreeE impact parameter (b),
105 // # of spectators n and p (fSpecn, fSpecp)
108 // // Get AliRun object from file or create it if not on file
109 // gAlice = (AliRun*)fBgrFile->Get("gAlice");
111 // gAlice = (AliRun*)fBgrFile->Get("gAlice");
112 // if (gAlice) printf("AliRun object found on file\n");
114 // printf("\n create new gAlice object");
115 // gAlice = new AliRun("gAlice","Alice test program");
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();
126 printf("\n HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
127 fNEvBgr,fImpPar,fSpecn,fSpecp);
130 //____________________________________________________________________________
131 TFile* AliZDCMerger::OpenBgrFile()
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);
143 //____________________________________________________________________________
144 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
145 Int_t &fFreeSpn, Int_t &fFreeSpp)
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++){
155 // Fragments generation
156 frag->GenerateIMF(zz, nAlpha);
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 2*nAlpha = %d, ztot = %d, ntot = %d\n",2*nAlpha, ztot, ntot);
167 printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
170 //____________________________________________________________________________
171 void AliZDCMerger::Mixing()
174 //printf("\n AliZDCMerger->Mixing\n");
176 // ### Background event Hits ###########################################
180 AliZDC *zdc = (AliZDC *)gAlice->GetModule("ZDC");
181 // if(zdc) printf("\n Ho trovato lo ZDC!\n");
183 // fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
185 char treeBgrName[20];
186 sprintf(treeBgrName,"TreeH%d",fNEvBgr);
187 fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
189 printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
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);
203 Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
204 //printf("\n --- ntracks = %d\n\n", ntracks);
206 Int_t itrack, nhits, ihit, j, sector[2];
208 AliZDCMergedHit *mHit;
210 TClonesArray &sdigits = *fMHits; // SDigits TCArray
214 for(itrack=0; itrack<ntracks; itrack++){
215 fTrHBgr->GetEvent(itrack);
217 nhits = fHitsBgr->GetEntries();
218 // printf(" nhits = %d \n", nhits);
219 for(ihit=0; ihit<nhits; ihit++){
220 zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
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);
232 new((*fMHits)[fNMhits]) AliZDCMergedHit(*mHit);
233 new (sdigits[fNMhits]) AliZDCMergedHit(*mHit);
239 //printf(" fNMhits (after bckg) = %d, \n",fNMhits);
241 // ### Signal event Hits ###########################################
247 //printf(" fNMhits (after signal) = %d \n",fNMhits);
251 //____________________________________________________________________________
252 void AliZDCMerger::ExtractSignal(Int_t SpecType)
255 // printf("\n Entering in Extract Signal\n");
258 if(SpecType == 1){ // --- Signal for spectator neutrons
259 fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
260 fSpecnFile = TFile::Open(fFnSpecn,"R");
262 printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn);
263 numEvents = fFreeSpn;
265 else if(SpecType == 2){ // --- Signal for spectator protons
266 fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
267 fSpecpFile = TFile::Open(fFnSpecp,"R");
269 printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp);
270 numEvents = fFreeSpp;
272 //printf("\n # of free spectators = %d\n", numEvents);
273 //printf("\n fNMhits (before adding signal) = %d\n",fNMhits);
275 TNtuple *zdcSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
276 Int_t nentries = (Int_t) zdcSignal->GetEntries();
277 //printf("\n # entries = %d\n", nentries);
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++){
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]);
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]);
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];
302 hitsSpec[j] = entry[j+3];
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);
309 else if(entry[0] > rnd[iev]){
314 }while(iev<numEvents);
317 //printf("\n fNMhits (after n signal) = %d\n",fNMhits);
320 else if(SpecType == 2){
321 //printf("\n fNMhits (after p signal) = %d\n",fNMhits);