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