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