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