1 /**************************************************************************
2 * Copyright(c) 1998-2000, 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 **************************************************************************/
18 #include <Riostream.h>
21 #include <TObjArray.h>
23 #include <TDirectory.h>
24 #include <TParticle.h>
26 #include "AliRICHMerger.h"
27 #include "AliRICHChamber.h"
28 #include "AliHitMap.h"
29 #include "AliRICHHitMapA1.h"
31 #include "AliRICHSDigit.h"
32 #include "AliRICHDigit.h"
33 #include "AliRICHTransientDigit.h"
37 ClassImp(AliRICHMerger)
39 //___________________________________________
40 AliRICHMerger::AliRICHMerger()
42 // Default constructor
57 //------------------------------------------------------------------------
58 AliRICHMerger::~AliRICHMerger()
61 if (fTrH1) delete fTrH1;
62 if (fHitsBgr) delete fHitsBgr;
63 if (fSDigitsBgr) delete fSDigitsBgr;
64 if (fHitMap) delete fHitMap;
65 if (fList) delete fList;
66 if (fBgrFile) delete fBgrFile;
69 //------------------------------------------------------------------------
70 Bool_t AliRICHMerger::Exists(const AliRICHSDigit *padhit)
72 return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
75 //------------------------------------------------------------------------
76 void AliRICHMerger::Update(AliRICHSDigit *padhit)
78 AliRICHTransientDigit *pdigit =
79 static_cast<AliRICHTransientDigit*>(
80 fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
84 Int_t iqpad = Int_t(padhit->QPad()); // charge per pad
85 pdigit->AddSignal(iqpad);
86 pdigit->AddPhysicsSignal(iqpad);
88 // update list of tracks
98 pdigit->UpdateTrackList(track,charge);
101 //------------------------------------------------------------------------
102 void AliRICHMerger::CreateNew(AliRICHSDigit *padhit)
104 fList->AddAtAndExpand(
105 new AliRICHTransientDigit(fNch,fDigits),fCounter);
106 fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
108 AliRICHTransientDigit* pdigit =
109 static_cast<AliRICHTransientDigit*>(fList->Last());
114 charge = padhit->QPad();
119 pdigit->AddToTrackList(track,charge);
124 //------------------------------------------------------------------------
125 void AliRICHMerger::Init()
129 if (fMerge && !fBgrFile) fBgrFile = InitBgr();
134 //------------------------------------------------------------------------
135 TFile* AliRICHMerger::InitBgr()
137 // Initialise background event
138 TFile *file = new TFile(fFnBgr);
139 // add error checking later
140 printf("\n AliRICHMerger has opened %s file with background event \n", fFnBgr);
141 fHitsBgr = new TClonesArray("AliRICHhit",1000);
142 fSDigitsBgr = new TClonesArray("AliRICHSDigit",1000);
146 //------------------------------------------------------------------------
147 void AliRICHMerger::Digitise(Int_t nev, Int_t flag)
150 // keep galice.root for signal and name differently the file for
151 // background when add! otherwise the track info for signal will be lost !
155 AliRICHChamber* iChamber;
156 AliSegmentation* segmentation;
158 fList = new TObjArray;
160 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
161 fHitMap= new AliHitMap* [kNCH];
166 // Get Hits Tree header from file
167 if(fHitsBgr) fHitsBgr->Clear();
168 if(fSDigitsBgr) fSDigitsBgr->Clear();
169 if(fTrH1) delete fTrH1;
173 sprintf(treeName,"TreeH%d",fEvNrBgr);
174 fTrH1 = (TTree*)gDirectory->Get(treeName);
176 printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
179 // Set branch addresses
182 sprintf(branchname,"%s",pRICH->GetName());
183 if (fTrH1 && fHitsBgr) {
184 branch = fTrH1->GetBranch(branchname);
185 if (branch) branch->SetAddress(&fHitsBgr);
187 if (fTrH1 && fSDigitsBgr) {
188 branch = fTrH1->GetBranch("RICHSDigits");
189 if (branch) branch->SetAddress(&fSDigitsBgr);
193 for (Int_t i =0; i<kNCH; i++) {
194 iChamber= &(pRICH->Chamber(i));
195 segmentation=iChamber->GetSegmentationModel(1);
196 fHitMap[i] = new AliRICHHitMapA1(segmentation, fList);
204 TTree *treeH = pRICH->TreeH();
205 Int_t ntracks =(Int_t) treeH->GetEntries();
206 for (fTrack=0; fTrack<ntracks; fTrack++) {
208 treeH->GetEvent(fTrack);
211 for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1);
213 mHit=(AliRICHhit*)pRICH->NextHit())
216 fNch = mHit->Chamber()-1; // chamber number
217 Int_t trackIndex = mHit->Track();
219 cerr<<"AliRICHMerger: chamber nr. fNch out of range: "<<fNch<<endl;
220 cerr<<" track: "<<fTrack<<endl;
223 iChamber = &(pRICH->Chamber(fNch));
227 // If flag is set, create digits only for some particles
229 TParticle *current = (TParticle*)gAlice->Particle(trackIndex);
230 if (current->GetPdgCode() >= 50000050)
232 TParticle *motherofcurrent = (TParticle*)gAlice->Particle(current->GetFirstMother());
233 particle = motherofcurrent->GetPdgCode();
237 particle = current->GetPdgCode();
241 //printf("Flag:%d\n",flag);
242 //printf("Track:%d\n",track);
243 //printf("Particle:%d\n",particle);
248 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
252 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
253 || TMath::Abs(particle)==311)
256 if (flag == 3 && TMath::Abs(particle)==2212)
259 if (flag == 4 && TMath::Abs(particle)==13)
262 if (flag == 5 && TMath::Abs(particle)==11)
265 if (flag == 6 && TMath::Abs(particle)==2112)
269 //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise);
276 // Loop over pad hits
277 for (AliRICHSDigit* mPad=
278 (AliRICHSDigit*)pRICH->FirstPad(mHit,pRICH->SDigits());
280 mPad=(AliRICHSDigit*)pRICH->NextPad(pRICH->SDigits()))
282 Int_t iqpad = mPad->QPad(); // charge per pad
283 fDigits[0]=mPad->PadX();
284 fDigits[1]=mPad->PadY();
287 fDigits[4]=mPad->HitNumber();
289 // build the list of fired pads and update the info
295 } //end loop over clusters
296 }// track type condition
300 // open the file with background
304 ntracks = (Int_t)fTrH1->GetEntries();
308 for (fTrack = 0; fTrack < ntracks; fTrack++) {
310 if (fHitsBgr) fHitsBgr->Clear();
311 if (fSDigitsBgr) fSDigitsBgr->Clear();
313 fTrH1->GetEvent(fTrack);
317 for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i)
319 mHit = (AliRICHhit*) (*fHitsBgr)[i];
320 fNch = mHit->Chamber()-1; // chamber number
321 iChamber = &(pRICH->Chamber(fNch));
324 // Loop over pad hits
325 for (AliRICHSDigit* mPad =
326 (AliRICHSDigit*)pRICH->FirstPad(mHit,fSDigitsBgr);
328 mPad = (AliRICHSDigit*)pRICH->NextPad(fSDigitsBgr))
330 fDigits[0] = mPad->PadX();
331 fDigits[1] = mPad->PadY();
332 fDigits[2] = mPad->QPad();
334 fDigits[4] = -1; // set hit number to -1 for bgr
336 // build the list of fired pads and update the info
342 } //end loop over clusters
346 TTree *treeK = gAlice->TreeK();
349 if (treeK) file = treeK->GetCurrentFile();
353 Int_t tracks[kMAXTRACKSPERRICHDIGIT];
354 Int_t charges[kMAXTRACKSPERRICHDIGIT];
355 Int_t nentries=fList->GetEntriesFast();
357 for (Int_t nent=0;nent<nentries;nent++) {
358 AliRICHTransientDigit *transDigit=(AliRICHTransientDigit*)fList->At(nent);
359 if (transDigit==0) continue;
360 Int_t ich=transDigit->GetChamber();
361 Int_t q=transDigit->Signal();
362 iChamber=&(pRICH->Chamber(ich));
363 AliRICHResponse * response=iChamber->GetResponseModel();
364 Int_t adcmax= (Int_t) response->MaxAdc();
367 // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
368 //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
369 Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());
371 //printf("Pedestal:%d\n",pedestal);
373 Float_t treshold = (pedestal + 4*2.2);
375 Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
376 Float_t noise = gRandom->Gaus(0, meanNoise);
377 q+=(Int_t)(noise + pedestal);
379 // magic number to be parametrised !!!
386 if ( q >= adcmax) q=adcmax;
387 fDigits[0]=transDigit->PadX();
388 fDigits[1]=transDigit->PadY();
390 fDigits[3]=transDigit->Physics();
391 fDigits[4]=transDigit->Hit();
393 Int_t nptracks = transDigit->GetNTracks();
395 // this was changed to accomodate the real number of tracks
396 if (nptracks > kMAXTRACKSPERRICHDIGIT) {
397 printf("Attention - tracks > 10 %d\n",nptracks);
398 nptracks=kMAXTRACKSPERRICHDIGIT;
401 printf("Attention - tracks > 2 %d \n",nptracks);
402 //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
403 //icat,ich,digits[0],digits[1],q);
405 for (Int_t tr=0;tr<nptracks;tr++) {
406 tracks[tr]=transDigit->GetTrack(tr);
407 charges[tr]=transDigit->GetCharge(tr);
408 } //end loop over list of tracks for one pad
409 if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
410 for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
417 //fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
420 pRICH->AddDigits(ich,tracks,charges,fDigits);
422 gAlice->TreeD()->Fill();
425 for(Int_t ii=0;ii<kNCH;++ii) {
432 TClonesArray *richDigits;
433 for (Int_t k=0;k<kNCH;k++) {
434 richDigits = pRICH->DigitsAddress(k);
435 Int_t ndigit=richDigits->GetEntriesFast();
436 printf ("Chamber %d digits %d \n",k,ndigit);
438 pRICH->ResetDigits(); /// ??? should it be here???
439 gAlice->TreeD()->Write(0,TObject::kOverwrite);
441 // gAlice->TreeD()->Reset();
442 // delete fList; // deleted in dtor
443 // gObjectTable->Print();