]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - MUON/AliMUONDigitizerv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Code for simulation, sdigitization and digitization moved from macros to compiled...
[u/mrichter/AliRoot.git] / MUON / AliMUONDigitizerv1.cxx
... / ...
CommitLineData
1
2
3#include <Riostream.h>
4#include <TDirectory.h>
5#include <TFile.h>
6#include <TObjArray.h>
7#include <TPDGCode.h>
8#include <TTree.h>
9#include <TMath.h>
10
11#include "AliMUON.h"
12#include "AliMUONChamber.h"
13#include "AliMUONConstants.h"
14#include "AliMUONDigit.h"
15#include "AliMUONDigitizerv1.h"
16#include "AliMUONHit.h"
17#include "AliMUONHitMapA1.h"
18#include "AliMUONPadHit.h"
19#include "AliMUONTransientDigit.h"
20#include "AliRun.h"
21#include "AliRunDigitizer.h"
22#include "AliRunLoader.h"
23#include "AliLoader.h"
24
25ClassImp(AliMUONDigitizerv1)
26
27//___________________________________________
28AliMUONDigitizerv1::AliMUONDigitizerv1() :
29 AliDigitizer(),
30 fHitMap(0),
31 fTDList(0),
32 fTDCounter(0),
33 fDebug(0),
34 fMask(0),
35 fSignal(0)
36{
37// Default ctor - don't use it
38 if (GetDebug()>2)
39 cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
40 <<"(AliRunDigitizer* manager) was processed"<<endl;
41}
42
43//___________________________________________
44AliMUONDigitizerv1::AliMUONDigitizerv1(AliRunDigitizer* manager):
45 AliDigitizer(manager),
46 fHitMap(0),
47 fTDList(0),
48 fTDCounter(0),
49 fDebug(0),
50 fMask(0),
51 fSignal(0)
52{
53// ctor which should be used
54 if (GetDebug()>2)
55 cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
56 <<"(AliRunDigitizer* manager) was processed"<<endl;
57}
58
59//------------------------------------------------------------------------
60AliMUONDigitizerv1::~AliMUONDigitizerv1()
61{
62// Destructor
63}
64
65//------------------------------------------------------------------------
66void AliMUONDigitizerv1::AddTransientDigit(AliMUONTransientDigit * mTD)
67{
68 // Choosing the maping of the cathode plane of the chamber:
69 Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
70 fTDList->AddAtAndExpand(mTD, fTDCounter);
71 fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
72 fTDCounter++;
73}
74
75//------------------------------------------------------------------------
76Bool_t AliMUONDigitizerv1::ExistTransientDigit(AliMUONTransientDigit * mTD)
77{
78 // Choosing the maping of the cathode plane of the chamber:
79 Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
80 return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
81}
82
83//------------------------------------------------------------------------
84Bool_t AliMUONDigitizerv1::Init()
85{
86
87// Initialization
88 if (GetDebug()>2) Info("Init","AliMUONDigitizerv1::Init() starts");
89
90 //Loaders (We assume input0 to be the output too)
91 AliRunLoader * runloader; // Input loader
92 AliLoader * gime;
93
94 // Getting runloader
95 runloader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
96 if (runloader == 0x0) {
97 Error("Init","RunLoader is not in input file 0");
98 return kFALSE; // RunDigitizer is not working.
99 }
100 // Getting MUONloader
101 gime = runloader->GetLoader("MUONLoader");
102 gime->LoadHits("READ");
103 gime->LoadDigits("RECREATE");
104
105 return kTRUE;
106}
107
108//------------------------------------------------------------------------
109void AliMUONDigitizerv1::UpdateTransientDigit(Int_t track, AliMUONTransientDigit * mTD)
110{
111 // Choosing the maping of the cathode plane of the chamber:
112 Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
113 AliMUONTransientDigit *pdigit =
114 static_cast<AliMUONTransientDigit*>(fHitMap[iNchCpl]->GetHit(mTD->PadX(),mTD->PadY()));
115 // update charge
116 //
117 Int_t iqpad = mTD->Signal(); // charge per pad
118 pdigit->AddSignal(iqpad);
119 pdigit->AddPhysicsSignal(iqpad);
120 // update list of tracks
121 //
122 Int_t charge;
123 track=+ fMask;
124 if (fSignal) charge = iqpad;
125 //else charge = kBgTag;
126 else charge = iqpad + fMask;
127
128 pdigit->UpdateTrackList(track,charge);
129}
130
131
132//--------------------------------------------------------------------------
133void AliMUONDigitizerv1::MakeTransientDigit(Int_t track, Int_t iHit, AliMUONHit * mHit)
134{
135 AliMUON *pMUON = (AliMUON *) gAlice->GetModule("MUON");
136 if (!pMUON) {
137 cerr<<"AliMUONDigitizerv1::Digitize Error:"
138 <<" module MUON not found in the input file"<<endl;
139 }
140 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit starts"<<endl;
141 Int_t ichamber = mHit->Chamber()-1;
142 AliMUONChamber & chamber = pMUON->Chamber(ichamber);
143 Float_t xhit = mHit->X();
144 Float_t yhit = mHit->Y();
145 Float_t zhit = mHit->Z();
146 Float_t eloss= mHit->Eloss();
147 Float_t tof = mHit->Age();
148 // Variables for chamber response from AliMUONChamber::DisIntegration
149 Float_t newdigit[6][500]; // Pad information
150 Int_t nnew=0; // Number of touched Pads per hit
151 Int_t digits[6];
152
153 //
154 // Calls the charge disintegration method of the current chamber
155 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit calling AliMUONChamber::DisIngtegration starts"<<endl;
156 chamber.DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newdigit);
157 // Creating a new TransientDigits from hit
158 for(Int_t iTD=0; iTD<nnew; iTD++) {
159 digits[0] = Int_t(newdigit[1][iTD]); // Padx of the Digit
160 digits[1] = Int_t(newdigit[2][iTD]); // Pady of the Digit
161 digits[2] = Int_t(newdigit[5][iTD]); // Cathode plane
162 digits[3] = Int_t(newdigit[3][iTD]); // Induced charge in the Pad
163 if (fSignal) digits[4] = Int_t(newdigit[3][iTD]);
164 else digits[4] = 0;
165 digits[5] = iHit+fMask; // Hit number in the list
166 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit " <<
167 "PadX "<< digits[0] << " " <<
168 "PadY "<< digits[1] << " " <<
169 "Plane " << digits[2] << " " <<
170 "Charge " << digits[3] <<" " <<
171 "Hit " << digits[5] << endl;
172 // list of tracks
173 Int_t charge;
174 track += fMask;
175 if (fSignal) charge = digits[3];
176 //else charge = kBgTag;
177 else charge = digits[3] + fMask;
178
179 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Creating AliMUONTransientDigit"<<endl;
180 AliMUONTransientDigit * mTD = new AliMUONTransientDigit(ichamber, digits);
181 mTD->AddToTrackList(track,charge);
182 if (!ExistTransientDigit(mTD)) {
183 AddTransientDigit(mTD);
184 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Adding TransientDigit"<<endl;
185 }
186 else {
187 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit updating TransientDigit"<<endl;
188 UpdateTransientDigit(track, mTD);
189 delete mTD;
190 }
191 }
192}
193//-----------------------------------------------------------------------
194void AliMUONDigitizerv1::Exec(Option_t* option)
195{
196 TString optionString = option;
197 if (optionString.Data() == "deb") {
198 Info("Digitize","Called with option deb ");
199 fDebug = 3;
200 }
201
202 AliMUONChamber* chamber;
203 AliSegmentation* c1Segmentation; //Cathode plane c1 of the chamber
204 AliSegmentation* c2Segmentation; //Cathode place c2 of the chamber
205
206 if (GetDebug()>2) Info("Digitize","AliMUONDigitizerv1::Digitize() starts");
207 fTDList = new TObjArray;
208
209 //Loaders (We assume input0 to be the output too)
210 AliRunLoader * runloader; // Input loader
211 AliLoader * gime;
212
213 // Getting runloader
214 runloader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
215 if (runloader == 0x0) {
216 Error("Digitize","RunLoader is not in input file 0");
217 return; // RunDigitizer is not working.
218 }
219 // Getting MUONloader
220 gime = runloader->GetLoader("MUONLoader");
221 if (gime->TreeH()==0x0) {
222 if (GetDebug()>2) Info("Digitize","TreeH is not loaded yet. Loading...");
223 gime->LoadHits("READ");
224 if (GetDebug()>2) Info("Digitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
225 }
226
227 if (GetDebug()>2) Info("Digitize","Loaders ready");
228
229 if (runloader->GetAliRun() == 0x0) runloader->LoadgAlice();
230 gAlice = runloader->GetAliRun();
231
232 // Getting Module MUON
233 AliMUON *pMUON = (AliMUON *) gAlice->GetDetector("MUON");
234 if (!pMUON) {
235 Error("Digitize","Module MUON not found in the input file");
236 return;
237 }
238 // Getting Muon data
239 AliMUONData * muondata = pMUON->GetMUONData();
240 muondata->SetLoader(gime);
241 muondata->SetTreeAddress("H");
242
243 // Loading Event
244 Int_t currentevent = fManager->GetOutputEventNr();
245
246 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() Event Number is "<<currentevent <<endl;
247 if ( (currentevent<10) ||
248 (Int_t(TMath::Log10(currentevent)) == TMath::Log10(currentevent) ) )
249 cout <<"ALiMUONDigitizerv1::Digitize() Event Number is "<< currentevent <<endl;
250
251 // Output file for digits same as hits
252 // AliRunLoader * runloaderout = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
253 //AliLoader * gimeout = runloaderout->GetLoader("MUONLoader");
254 // New branch per chamber for MUON digit in the tree of digits
255 if (gime->TreeD() == 0x0) {
256 gime->MakeDigitsContainer();
257 }
258 TTree* treeD = gime->TreeD();
259 muondata->MakeBranch("D");
260 muondata->SetTreeAddress("D");
261
262 // Array of pointer of the AliMUONHitMapA1:
263 // two HitMaps per chamber, or one HitMap per cahtode plane
264 fHitMap= new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];
265
266 //Loop over chambers for the definition AliMUONHitMap
267 for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
268 chamber = &(pMUON->Chamber(i));
269 c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
270 fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
271 c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
272 fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
273 }
274
275// Loop over files to merge and to digitize
276 fSignal = kTRUE;
277 for (Int_t inputFile=0; inputFile<fManager->GetNinputs(); inputFile++) {
278 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() Input File is "<<inputFile<<endl;
279
280
281 // Connect MUON Hit branch
282 if (inputFile > 0 ) {
283 fSignal = kFALSE;
284 runloader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
285 if (runloader == 0x0) {
286 cerr<<"AliMUONDigitizerv1::Digitize() RunLoader for inputFile "<<inputFile<< " not found !!! "<<endl;
287 }
288 gime = runloader->GetLoader("MUONLoader");
289 if (gime->TreeH() == 0x0) gime->LoadHits("READ");
290 muondata->SetLoader(gime);
291 muondata->SetTreeAddress("H");
292 }
293
294 // Setting the address
295 TTree *treeH = gime->TreeH();
296 if (treeH == 0x0) {
297 Error("Digitize","Can not get TreeH from input %d",inputFile);
298 Info("Digitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
299 return;
300 }
301 if (GetDebug()>2) {
302 cerr<<"AliMUONDigitizerv1::Exec inputFile is "<<inputFile<<" "<<endl;
303 cerr<<"AliMUONDigitizerv1::Exec treeH" << treeH <<endl;
304 }
305
306 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec Setting tree addresses"<<endl;
307
308 fMask = fManager->GetMask(inputFile);
309 //
310 // Loop over tracks
311 Int_t itrack;
312 Int_t ntracks = (Int_t) treeH->GetEntries();
313 for (itrack = 0; itrack < ntracks; itrack++) {
314 if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec itrack = "<<itrack<<endl;
315 muondata->ResetHits();
316 treeH->GetEvent(itrack);
317 //
318 // Loop over hits
319 Int_t ihit, ichamber;
320 AliMUONHit* mHit;
321 TClonesArray* hits = muondata->Hits();
322 for(ihit = 0; ihit < hits->GetEntriesFast(); ihit++) {
323 mHit = static_cast<AliMUONHit*>(hits->At(ihit));
324 ichamber = mHit->Chamber()-1; // chamber number
325 if (ichamber > AliMUONConstants::NCh()-1) {
326 cerr<<"AliMUONDigitizer: ERROR: "
327 <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
328 <<ichamber<<", "<< AliMUONConstants::NCh()<<endl;
329 return;
330 }
331 chamber = &(pMUON->Chamber(ichamber));
332 //
333 //Dumping Hit content:
334 if (GetDebug()>2) {
335 cerr<<"AliMuonDigitizerv1::Exec ihit, ichamber, x, y, z, eloss " <<
336 ihit << " " <<
337 mHit->Chamber() << " " <<
338 mHit->X() << " " <<
339 mHit->Y() << " " <<
340 mHit->Z() << " " <<
341 mHit->Eloss() << " " << endl;
342 }
343 //
344 // Inititializing Correlation
345 chamber->ChargeCorrelationInit();
346 if (ichamber < AliMUONConstants::NTrackingCh()) {
347 // Tracking Chamber
348 // Initialize hit position (cursor) in the segmentation model
349 chamber->SigGenInit(mHit->X(), mHit->Y(), mHit->Z());
350 } else {
351 // Trigger Chamber
352 }
353 MakeTransientDigit(itrack, ihit, mHit);
354 } // hit loop
355 } // track loop
356 } // end file loop
357 if (GetDebug()>2) cerr<<"AliMUONDigitizer::Exec End of hits, track and file loops"<<endl;
358
359 // Loop on cathodes
360 Int_t icat;
361 for(icat=0; icat<2; icat++) {
362 //
363 // Filling Digit List
364 Int_t tracks[kMAXTRACKS];
365 Int_t charges[kMAXTRACKS];
366 Int_t nentries = fTDList->GetEntriesFast();
367 Int_t digits[6];
368 for (Int_t nent = 0; nent < nentries; nent++) {
369 AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
370 if (address == 0) continue;
371 Int_t ich = address->Chamber();
372 Int_t q = address->Signal();
373 chamber = &(pMUON->Chamber(ich));
374 //
375 // Digit Response (noise, threshold, saturation, ...)
376 AliMUONResponse * response = chamber->ResponseModel();
377 q = response->DigitResponse(q,address);
378
379 if (!q) continue;
380
381 digits[0] = address->PadX();
382 digits[1] = address->PadY();
383 digits[2] = address->Cathode()-1;
384 digits[3] = q;
385 digits[4] = address->Physics();
386 digits[5] = address->Hit();
387
388 Int_t nptracks = address->GetNTracks();
389
390 if (nptracks > kMAXTRACKS) {
391 if (GetDebug() >0) {
392 cerr<<"AliMUONDigitizer:Exec nptracks > 10 "<<nptracks;
393 cerr<<"reset to max value "<<kMAXTRACKS<<endl;
394 }
395 nptracks = kMAXTRACKS;
396 }
397 if (nptracks > 2 && GetDebug() >2) {
398 cerr<<"AliMUONDigitizer::Exec nptracks > 2 "<<nptracks<<endl;
399 // printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
400 }
401 for (Int_t tr = 0; tr < nptracks; tr++) {
402 tracks[tr] = address->GetTrack(tr);
403 charges[tr] = address->GetCharge(tr);
404 } //end loop over list of tracks for one pad
405 // Sort list of tracks according to charge
406 if (nptracks > 1) {
407 SortTracks(tracks,charges,nptracks);
408 }
409 if (nptracks < kMAXTRACKS ) {
410 for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
411 tracks[i] = -1;
412 charges[i] = 0;
413 }
414 }
415
416 // Add digits
417 if (GetDebug()>3) cerr<<"AliMUONDigitzerv1::Exex TransientDigit to Digit"<<endl;
418 if ( digits[2] == icat ) muondata->AddDigit(ich,tracks,charges,digits);
419// printf("test rm ich %d padX %d padY %d \n",ich, digits[0], digits[1]);
420 }
421 // Filling list of digits per chamber for a given cathode.
422 muondata->Fill("D");
423 muondata->ResetDigits();
424 } // end loop cathode
425 fTDList->Delete();
426
427 for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
428 if (fHitMap[ii]) {
429 delete fHitMap[ii];
430 fHitMap[ii] = 0;
431 }
432 }
433
434 if (GetDebug()>2)
435 cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
436 <<treeD->GetName()<<endl;
437
438 gime->WriteDigits("OVERWRITE");
439 delete [] fHitMap;
440 delete fTDList;
441 muondata->ResetHits();
442 gime->UnloadHits();
443 gime->UnloadDigits();
444}
445//------------------------------------------------------------------------
446void AliMUONDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
447{
448 //
449 // Sort the list of tracks contributing to a given digit
450 // Only the 3 most significant tracks are acctually sorted
451 //
452
453 //
454 // Loop over signals, only 3 times
455 //
456
457 Int_t qmax;
458 Int_t jmax;
459 Int_t idx[3] = {-2,-2,-2};
460 Int_t jch[3] = {-2,-2,-2};
461 Int_t jtr[3] = {-2,-2,-2};
462 Int_t i,j,imax;
463
464 if (ntr<3) imax=ntr;
465 else imax=3;
466 for(i=0;i<imax;i++){
467 qmax=0;
468 jmax=0;
469
470 for(j=0;j<ntr;j++){
471
472 if((i == 1 && j == idx[i-1])
473 ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
474
475 if(charges[j] > qmax) {
476 qmax = charges[j];
477 jmax=j;
478 }
479 }
480
481 if(qmax > 0) {
482 idx[i]=jmax;
483 jch[i]=charges[jmax];
484 jtr[i]=tracks[jmax];
485 }
486
487 }
488
489 for(i=0;i<3;i++){
490 if (jtr[i] == -2) {
491 charges[i]=0;
492 tracks[i]=0;
493 } else {
494 charges[i]=jch[i];
495 tracks[i]=jtr[i];
496 }
497 }
498}