]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RICH/AliRICHMerger.cxx
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / RICH / AliRICHMerger.cxx
CommitLineData
d8a72780 1/**************************************************************************
2 * Copyright(c) 1998-2000, 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$Log$
adf7bb67 18Revision 1.7 2001/11/02 15:37:26 hristov
19Digitizer class created. Code cleaning and bug fixes (J.Chudoba)
20
31eeed4f 21Revision 1.5 2001/10/23 13:03:35 hristov
22The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba)
23
6e585aa2 24Revision 1.4 2001/10/21 18:31:24 hristov
25Several pointers were set to zero in the default constructors to avoid memory management problems
26
2685bf00 27Revision 1.3 2001/05/16 14:57:20 alibrary
28New files for folders and Stack
29
9e1a0ddb 30Revision 1.2 2001/03/14 18:16:08 jbarbosa
31Corrected bug (more to correct).
32File "points.dat" is no longer created.
33
14886cb7 34Revision 1.1 2001/02/27 22:13:34 jbarbosa
35Implementing merger class.
36
d8a72780 37*/
adf7bb67 38#include <iostream.h>
d8a72780 39
40#include <TTree.h>
d8a72780 41#include <TObjArray.h>
42#include <TFile.h>
43#include <TDirectory.h>
44#include <TParticle.h>
45
d8a72780 46#include "AliRICHMerger.h"
47#include "AliRICHChamber.h"
48#include "AliHitMap.h"
49#include "AliRICHHitMapA1.h"
50#include "AliRICH.h"
51#include "AliRICHHit.h"
52#include "AliRICHSDigit.h"
53#include "AliRICHDigit.h"
54#include "AliRICHTransientDigit.h"
55#include "AliRun.h"
56#include "AliPDG.h"
57
58ClassImp(AliRICHMerger)
59
60//___________________________________________
b762c2f6 61 AliRICHMerger::AliRICHMerger()
d8a72780 62{
63// Default constructor
b762c2f6 64 fEvNrSig = 0;
65 fEvNrBgr = 0;
66 fMerge = kDigitize;
67 fFnBgr = 0;
68 fHitMap = 0;
69 fList = 0;
70 fTrH1 = 0;
71 fHitsBgr = 0;
72 fSDigitsBgr = 0;
73 fHitMap = 0;
74 fList = 0;
75 fBgrFile = 0;
d8a72780 76}
77
78//------------------------------------------------------------------------
79AliRICHMerger::~AliRICHMerger()
80{
81// Destructor
b762c2f6 82 if (fTrH1) delete fTrH1;
83 if (fHitsBgr) delete fHitsBgr;
84 if (fSDigitsBgr) delete fSDigitsBgr;
85 if (fHitMap) delete fHitMap;
86 if (fList) delete fList;
87 if (fBgrFile) delete fBgrFile;
d8a72780 88}
89
90//------------------------------------------------------------------------
b762c2f6 91Bool_t AliRICHMerger::Exists(const AliRICHSDigit *padhit)
d8a72780 92{
b762c2f6 93 return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
d8a72780 94}
95
96//------------------------------------------------------------------------
b762c2f6 97void AliRICHMerger::Update(AliRICHSDigit *padhit)
d8a72780 98{
b762c2f6 99 AliRICHTransientDigit *pdigit =
100 static_cast<AliRICHTransientDigit*>(
101 fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
102
103 // update charge
104 //
105 Int_t iqpad = Int_t(padhit->QPad()); // charge per pad
106 pdigit->AddSignal(iqpad);
107 pdigit->AddPhysicsSignal(iqpad);
108
109 // update list of tracks
110 //
111 Int_t track, charge;
112 if (fSignal) {
113 track = fTrack;
114 charge = iqpad;
115 } else {
116 track = kBgTag;
117 charge = kBgTag;
118 }
119 pdigit->UpdateTrackList(track,charge);
d8a72780 120}
121
122//------------------------------------------------------------------------
b762c2f6 123void AliRICHMerger::CreateNew(AliRICHSDigit *padhit)
d8a72780 124{
b762c2f6 125 fList->AddAtAndExpand(
126 new AliRICHTransientDigit(fNch,fDigits),fCounter);
127 fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
128
129 AliRICHTransientDigit* pdigit =
130 static_cast<AliRICHTransientDigit*>(fList->Last());
131 // list of tracks
132 Int_t track, charge;
133 if (fSignal) {
134 track = fTrack;
135 charge = padhit->QPad();
136 } else {
137 track = kBgTag;
138 charge = kBgTag;
139 }
140 pdigit->AddToTrackList(track,charge);
141 fCounter++;
d8a72780 142}
143
144
145//------------------------------------------------------------------------
146void AliRICHMerger::Init()
147{
148// Initialisation
149
b762c2f6 150 if (fMerge && !fBgrFile) fBgrFile = InitBgr();
d8a72780 151}
152
153
154
155//------------------------------------------------------------------------
156TFile* AliRICHMerger::InitBgr()
157{
158// Initialise background event
b762c2f6 159 TFile *file = new TFile(fFnBgr);
d8a72780 160// add error checking later
b762c2f6 161 printf("\n AliRICHMerger has opened %s file with background event \n", fFnBgr);
162 fHitsBgr = new TClonesArray("AliRICHHit",1000);
163 fSDigitsBgr = new TClonesArray("AliRICHSDigit",1000);
164 return file;
d8a72780 165}
166
167//------------------------------------------------------------------------
168void AliRICHMerger::Digitise(Int_t nev, Int_t flag)
169{
170
b762c2f6 171 // keep galice.root for signal and name differently the file for
172 // background when add! otherwise the track info for signal will be lost !
d8a72780 173
174 Int_t particle;
175
d8a72780 176 AliRICHChamber* iChamber;
177 AliSegmentation* segmentation;
178
b762c2f6 179 fList = new TObjArray;
d8a72780 180
b762c2f6 181 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
182 fHitMap= new AliHitMap* [kNCH];
d8a72780 183
b762c2f6 184 if (fMerge ) {
185 fBgrFile->cd();
186 //
187 // Get Hits Tree header from file
188 if(fHitsBgr) fHitsBgr->Clear();
189 if(fSDigitsBgr) fSDigitsBgr->Clear();
190 if(fTrH1) delete fTrH1;
191 fTrH1 = 0;
d8a72780 192
b762c2f6 193 char treeName[20];
194 sprintf(treeName,"TreeH%d",fEvNrBgr);
195 fTrH1 = (TTree*)gDirectory->Get(treeName);
196 if (!fTrH1) {
197 printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
d8a72780 198 }
d8a72780 199 //
b762c2f6 200 // Set branch addresses
201 TBranch *branch;
202 char branchname[20];
203 sprintf(branchname,"%s",pRICH->GetName());
204 if (fTrH1 && fHitsBgr) {
205 branch = fTrH1->GetBranch(branchname);
206 if (branch) branch->SetAddress(&fHitsBgr);
207 }
208 if (fTrH1 && fSDigitsBgr) {
209 branch = fTrH1->GetBranch("RICHSDigits");
210 if (branch) branch->SetAddress(&fSDigitsBgr);
211 }
212 }
d8a72780 213
b762c2f6 214 for (Int_t i =0; i<kNCH; i++) {
215 iChamber= &(pRICH->Chamber(i));
216 segmentation=iChamber->GetSegmentationModel(1);
217 fHitMap[i] = new AliRICHHitMapA1(segmentation, fList);
218 }
219 //
220 // Loop over tracks
221 //
222
223 fSignal = kTRUE;
224 fCounter = 0;
225 TTree *treeH = gAlice->TreeH();
226 Int_t ntracks =(Int_t) treeH->GetEntries();
227 for (fTrack=0; fTrack<ntracks; fTrack++) {
228 gAlice->ResetHits();
229 treeH->GetEvent(fTrack);
230 //
231 // Loop over hits
232 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
233 mHit;
234 mHit=(AliRICHHit*)pRICH->NextHit())
235 {
d8a72780 236
b762c2f6 237 fNch = mHit->Chamber()-1; // chamber number
238 Int_t trackIndex = mHit->Track();
239 if (fNch >= kNCH) {
240 cerr<<"AliRICHMerger: chamber nr. fNch out of range: "<<fNch<<endl;
241 cerr<<" track: "<<fTrack<<endl;
242 continue;
243 }
244 iChamber = &(pRICH->Chamber(fNch));
d8a72780 245
b762c2f6 246
247//
248// If flag is set, create digits only for some particles
249//
250 TParticle *current = (TParticle*)gAlice->Particle(trackIndex);
251 if (current->GetPdgCode() >= 50000050)
252 {
253 TParticle *motherofcurrent = (TParticle*)gAlice->Particle(current->GetFirstMother());
254 particle = motherofcurrent->GetPdgCode();
255 }
256 else
257 {
258 particle = current->GetPdgCode();
259 }
d8a72780 260
261
b762c2f6 262 //printf("Flag:%d\n",flag);
263 //printf("Track:%d\n",track);
264 //printf("Particle:%d\n",particle);
d8a72780 265
b762c2f6 266 Int_t digitise=1;
d8a72780 267
b762c2f6 268 if (flag == 1)
269 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
270 digitise=0;
d8a72780 271
b762c2f6 272 if (flag == 2)
273 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
274 || TMath::Abs(particle)==311)
275 digitise=0;
d8a72780 276
b762c2f6 277 if (flag == 3 && TMath::Abs(particle)==2212)
278 digitise=0;
d8a72780 279
b762c2f6 280 if (flag == 4 && TMath::Abs(particle)==13)
281 digitise=0;
d8a72780 282
b762c2f6 283 if (flag == 5 && TMath::Abs(particle)==11)
284 digitise=0;
d8a72780 285
b762c2f6 286 if (flag == 6 && TMath::Abs(particle)==2112)
287 digitise=0;
d8a72780 288
289
b762c2f6 290 //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise);
d8a72780 291
292
b762c2f6 293 if (digitise)
294 {
d8a72780 295
b762c2f6 296 //
297 // Loop over pad hits
298 for (AliRICHSDigit* mPad=
299 (AliRICHSDigit*)pRICH->FirstPad(mHit,pRICH->SDigits());
300 mPad;
301 mPad=(AliRICHSDigit*)pRICH->NextPad(pRICH->SDigits()))
302 {
303 Int_t iqpad = mPad->QPad(); // charge per pad
304 fDigits[0]=mPad->PadX();
305 fDigits[1]=mPad->PadY();
306 fDigits[2]=iqpad;
307 fDigits[3]=iqpad;
308 fDigits[4]=mPad->HitNumber();
d8a72780 309
b762c2f6 310 // build the list of fired pads and update the info
311 if (Exists(mPad)) {
312 Update(mPad);
313 } else {
314 CreateNew(mPad);
315 }
316 } //end loop over clusters
317 }// track type condition
318 } // hit loop
319 } // track loop
d8a72780 320
b762c2f6 321 // open the file with background
d8a72780 322
b762c2f6 323 if (fMerge) {
324 fSignal = kFALSE;
325 ntracks = (Int_t)fTrH1->GetEntries();
326 //
327 // Loop over tracks
328 //
329 for (fTrack = 0; fTrack < ntracks; fTrack++) {
d8a72780 330
b762c2f6 331 if (fHitsBgr) fHitsBgr->Clear();
332 if (fSDigitsBgr) fSDigitsBgr->Clear();
d8a72780 333
b762c2f6 334 fTrH1->GetEvent(fTrack);
335 //
336 // Loop over hits
337 AliRICHHit* mHit;
338 for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
339 {
340 mHit = (AliRICHHit*) (*fHitsBgr)[i];
341 fNch = mHit->Chamber()-1; // chamber number
342 iChamber = &(pRICH->Chamber(fNch));
343
d8a72780 344 //
b762c2f6 345 // Loop over pad hits
346 for (AliRICHSDigit* mPad =
347 (AliRICHSDigit*)pRICH->FirstPad(mHit,fSDigitsBgr);
348 mPad;
349 mPad = (AliRICHSDigit*)pRICH->NextPad(fSDigitsBgr))
350 {
351 fDigits[0] = mPad->PadX();
352 fDigits[1] = mPad->PadY();
353 fDigits[2] = mPad->QPad();
354 fDigits[3] = 0;
355 fDigits[4] = -1; // set hit number to -1 for bgr
d8a72780 356
b762c2f6 357 // build the list of fired pads and update the info
358 if (!Exists(mPad)) {
359 CreateNew(mPad);
360 } else {
361 Update(mPad);
362 } // end if !Exists
363 } //end loop over clusters
364 } // hit loop
365 } // track loop
d8a72780 366
b762c2f6 367 TTree *treeK = gAlice->TreeK();
368 TFile *file = NULL;
d8a72780 369
b762c2f6 370 if (treeK) file = treeK->GetCurrentFile();
371 file->cd();
372 } // if fMerge
d8a72780 373
b762c2f6 374 Int_t tracks[kMAXTRACKSPERRICHDIGIT];
375 Int_t charges[kMAXTRACKSPERRICHDIGIT];
376 Int_t nentries=fList->GetEntriesFast();
d8a72780 377
b762c2f6 378 for (Int_t nent=0;nent<nentries;nent++) {
379 AliRICHTransientDigit *transDigit=(AliRICHTransientDigit*)fList->At(nent);
380 if (transDigit==0) continue;
381 Int_t ich=transDigit->GetChamber();
382 Int_t q=transDigit->Signal();
383 iChamber=&(pRICH->Chamber(ich));
384 AliRICHResponse * response=iChamber->GetResponseModel();
385 Int_t adcmax= (Int_t) response->MaxAdc();
d8a72780 386
d8a72780 387
b762c2f6 388 // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
389 //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
390 Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());
d8a72780 391
b762c2f6 392 //printf("Pedestal:%d\n",pedestal);
393 //Int_t pedestal=0;
394 Float_t treshold = (pedestal + 4*2.2);
d8a72780 395
b762c2f6 396 Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
397 Float_t noise = gRandom->Gaus(0, meanNoise);
398 q+=(Int_t)(noise + pedestal);
399 //q+=(Int_t)(noise);
400 // magic number to be parametrised !!!
401 if ( q <= treshold)
402 {
d8a72780 403 q = q - pedestal;
b762c2f6 404 continue;
405 }
406 q = q - pedestal;
407 if ( q >= adcmax) q=adcmax;
408 fDigits[0]=transDigit->PadX();
409 fDigits[1]=transDigit->PadY();
410 fDigits[2]=q;
411 fDigits[3]=transDigit->Physics();
412 fDigits[4]=transDigit->Hit();
413
414 Int_t nptracks = transDigit->GetNTracks();
415
416 // this was changed to accomodate the real number of tracks
417 if (nptracks > kMAXTRACKSPERRICHDIGIT) {
418 printf("Attention - tracks > 10 %d\n",nptracks);
419 nptracks=kMAXTRACKSPERRICHDIGIT;
420 }
421 if (nptracks > 2) {
422 printf("Attention - tracks > 2 %d \n",nptracks);
423 //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
424 //icat,ich,digits[0],digits[1],q);
425 }
426 for (Int_t tr=0;tr<nptracks;tr++) {
427 tracks[tr]=transDigit->GetTrack(tr);
428 charges[tr]=transDigit->GetCharge(tr);
429 } //end loop over list of tracks for one pad
430 if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
431 for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
432 tracks[t]=0;
433 charges[t]=0;
d8a72780 434 }
b762c2f6 435 }
436 //write file
437 //if (ich==2)
438 //fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
d8a72780 439
b762c2f6 440 // fill digits
441 pRICH->AddDigits(ich,tracks,charges,fDigits);
442 }
443 gAlice->TreeD()->Fill();
444
445 fList->Delete();
446 for(Int_t ii=0;ii<kNCH;++ii) {
447 if (fHitMap[ii]) {
448 delete fHitMap[ii];
449 fHitMap[ii]=0;
d8a72780 450 }
b762c2f6 451 }
d8a72780 452
b762c2f6 453 TClonesArray *richDigits;
454 for (Int_t k=0;k<kNCH;k++) {
455 richDigits = pRICH->DigitsAddress(k);
456 Int_t ndigit=richDigits->GetEntriesFast();
457 printf ("Chamber %d digits %d \n",k,ndigit);
458 }
459 pRICH->ResetDigits(); /// ??? should it be here???
460 gAlice->TreeD()->Write(0,TObject::kOverwrite);
461 // reset tree
462 // gAlice->TreeD()->Reset();
463 // delete fList; // deleted in dtor
464 // gObjectTable->Print();
d8a72780 465
466}