Code for simulation, sdigitization and digitization moved from macros to compiled...
[u/mrichter/AliRoot.git] / MUON / AliMUONDigitizerv1.cxx
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
25 ClassImp(AliMUONDigitizerv1)
26
27 //___________________________________________
28 AliMUONDigitizerv1::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 //___________________________________________
44 AliMUONDigitizerv1::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 //------------------------------------------------------------------------
60 AliMUONDigitizerv1::~AliMUONDigitizerv1()
61 {
62 // Destructor
63 }
64
65 //------------------------------------------------------------------------
66 void 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 //------------------------------------------------------------------------
76 Bool_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 //------------------------------------------------------------------------
84 Bool_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 //------------------------------------------------------------------------
109 void 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 //--------------------------------------------------------------------------
133 void 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 //-----------------------------------------------------------------------
194 void 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 //------------------------------------------------------------------------
446 void 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 }