95427b9dd31dd4afb164a43471ca78f75495c2a9
[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() :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 //___________________________________________
39 AliMUONDigitizerv1::AliMUONDigitizerv1(AliRunDigitizer* manager) 
40     :AliDigitizer(manager)
41 {
42 // ctor which should be used
43   fHitMap  = 0;
44   fTDList  = 0;
45   fDebug   = 0; 
46   fHits = new TClonesArray("AliMUONHit",1000);
47   if (GetDebug()>2) 
48     cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
49         <<"(AliRunDigitizer* manager) was processed"<<endl;
50 }
51
52 //------------------------------------------------------------------------
53 AliMUONDigitizerv1::~AliMUONDigitizerv1()
54 {
55 // Destructor
56   delete fHits;
57 }
58
59 //------------------------------------------------------------------------
60 void 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 //------------------------------------------------------------------------
70 Bool_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 //------------------------------------------------------------------------
78 Bool_t AliMUONDigitizerv1::Init()
79 {
80 // Initialization 
81   return kTRUE;
82 }
83
84 //------------------------------------------------------------------------
85 void 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   //
98   Int_t charge;   
99   track=+ fMask;
100   if (fSignal) charge = iqpad;
101   //else         charge = kBgTag;
102   else         charge = iqpad + fMask; 
103                
104   pdigit->UpdateTrackList(track,charge);
105 }
106
107
108 //--------------------------------------------------------------------------
109 void 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
149     Int_t charge;   
150     track += fMask;
151     if (fSignal) charge = digits[3];
152     //else         charge = kBgTag;
153     else         charge = digits[3] + fMask;
154                  
155     if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Creating AliMUONTransientDigit"<<endl;
156     AliMUONTransientDigit * mTD = new AliMUONTransientDigit(ichamber, digits);
157     mTD->AddToTrackList(track,charge);
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 //-----------------------------------------------------------------------
170 void AliMUONDigitizerv1::Exec(Option_t* option)
171 {
172   TString optionString = option;
173   if (optionString.Data() == "deb") {
174     cout<<"AliMUONDigitizerv1::Exec: called with option deb "<<endl;
175     fDebug = 3;
176   }
177   AliMUONChamber*   chamber;
178   AliSegmentation*  c1Segmentation; //Cathode plane c1 of the chamber
179   AliSegmentation*  c2Segmentation; //Cathode place c2 of the chamber
180
181   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() starts"<<endl;
182   fTDList = new TObjArray;
183
184   //Loaders (We assume input0 to be the output too)
185   AliRunLoader * runloader;    // Input   loader
186   AliLoader    * gime; 
187  
188   runloader    = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
189   if (runloader == 0x0) {
190     cerr<<"AliMUONDigitizerv1::Digitize() opening file "<<fManager->GetInputFileName(0,0)<<endl;
191     return;  // RunDigitizer is not working.  
192   }
193   gime         = runloader->GetLoader("MUONLoader");
194   if (gime->TreeH()==0x0) {
195      Info("Digitize","TreeH is not loaded yet. Loading...");
196      gime->LoadHits("READ");
197      Info("Digitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
198   }
199
200   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() loaders"<<endl;
201
202    if (runloader->GetAliRun() == 0x0) runloader->LoadgAlice();
203   gAlice = runloader->GetAliRun();
204
205   // Getting Module MUON  
206   AliMUON *pMUON  = (AliMUON *) gAlice->GetDetector("MUON");
207   if (!pMUON) {
208     cerr<<"AliMUONDigitizerv1::Digitize Error:"
209         <<" module MUON not found in the input file"<<endl;
210     return;
211   } 
212
213   // Loading Event
214   Int_t currentevent = fManager->GetOutputEventNr();
215   
216   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() Event Number is "<<currentevent <<endl;
217   if ( (currentevent<10)                                                 || 
218        (Int_t(TMath::Log10(currentevent)) == TMath::Log10(currentevent) ) )
219     cout <<"ALiMUONDigitizerv1::Digitize() Event Number is "<< currentevent <<endl;
220
221   // Output file for digits 
222   AliRunLoader *  runloaderout  = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
223   AliLoader * gimeout         = runloaderout->GetLoader("MUONLoader"); 
224   // New branch per chamber for MUON digit in the tree of digits
225   if (gime->TreeD() == 0x0) gimeout->MakeDigitsContainer();
226   TTree* treeD = gimeout->TreeD();
227   pMUON->MakeBranchInTreeD(treeD);
228
229   // Array of pointer of the AliMUONHitMapA1:
230   //  two HitMaps per chamber, or one HitMap per cahtode plane
231   fHitMap= new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];
232
233   //Loop over chambers for the definition AliMUONHitMap
234   for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
235     chamber = &(pMUON->Chamber(i));
236     c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
237     fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
238     c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
239     fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
240   }
241
242 // Loop over files to merge and to digitize
243     fSignal = kTRUE;
244     for (Int_t inputFile=0; inputFile<fManager->GetNinputs(); inputFile++) {
245       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() Input File is "<<inputFile<<endl;
246
247       // Hits Branch
248       TBranch *branchHits = 0;
249
250       // Connect MUON Hit branch
251       if (inputFile > 0 ) {
252         fSignal = kFALSE;
253          runloader    = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
254          if (runloader == 0x0) {
255            cerr<<"AliMUONDigitizerv1::Digitize() RunLoader for inputFile "<<inputFile<< " not found !!! "<<endl;
256          }
257          gime  = runloader->GetLoader("MUONLoader");
258          if (gime->TreeH() == 0x0) gime->LoadHits("READ");      
259       }
260
261       // Setting the address of fHits list
262       TTree *treeH = gime->TreeH();
263       if (treeH == 0x0) {
264         Error("Digitize","Can not get TreeH from input %d",inputFile);
265         Info("Digitize","Now treeH is %#x. MUONLoader is %#x",gime->TreeH(),gime);
266         return;
267       }
268       if (GetDebug()>2) {
269         cerr<<"AliMUONDigitizerv1::Exec inputFile is "<<inputFile<<" "<<endl;
270         cerr<<"AliMUONDigitizerv1::Exec treeH, fHits "<<treeH<<" "<<fHits<<endl;
271       }
272       if (treeH && fHits) {
273         branchHits = treeH->GetBranch("MUON");
274         if (branchHits) {
275           fHits->Clear();
276           branchHits->SetAddress(&fHits);
277         }
278         else
279           Error("Exec","branch MUON was not found");
280       }
281       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec branchHits = "<<branchHits<<endl;
282
283
284       fMask = fManager->GetMask(inputFile);
285       //
286       // Loop over tracks
287       Int_t itrack;
288       Int_t ntracks = (Int_t) treeH->GetEntries();
289       for (itrack = 0; itrack < ntracks; itrack++) {
290         if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec itrack = "<<itrack<<endl;
291         fHits->Clear();
292         treeH->GetEvent(itrack);
293         //
294         //  Loop over hits
295         Int_t ihit, ichamber;
296         AliMUONHit* mHit;
297         for(ihit = 0; ihit < fHits->GetEntriesFast(); ihit++) {
298           mHit = static_cast<AliMUONHit*>(fHits->At(ihit));
299           ichamber = mHit->Chamber()-1;  // chamber number
300           if (ichamber > AliMUONConstants::NCh()-1) {
301             cerr<<"AliMUONDigitizer: ERROR: "
302                 <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
303                 <<ichamber<<", "<< AliMUONConstants::NCh()<<endl;
304             return;
305           }
306           chamber = &(pMUON->Chamber(ichamber));
307           //
308           //Dumping Hit content:
309           if (GetDebug()>2) {
310             cerr<<"AliMuonDigitizerv1::Exec ihit, ichamber, x, y, z, eloss " <<
311               ihit << " " << 
312               mHit->Chamber() << " " <<
313               mHit->X() << " " <<
314               mHit->Y() << " " <<
315               mHit->Z() << " " <<
316               mHit->Eloss() << " " << endl;
317           }
318           // 
319           // Inititializing Correlation
320           chamber->ChargeCorrelationInit();
321           if (ichamber < AliMUONConstants::NTrackingCh()) {
322             // Tracking Chamber
323             // Initialize hit position (cursor) in the segmentation model 
324             chamber->SigGenInit(mHit->X(), mHit->Y(), mHit->Z());
325           } else {
326             // Trigger Chamber 
327           }
328           MakeTransientDigit(itrack, ihit, mHit);
329         } // hit loop
330       } // track loop
331     } // end file loop
332     if (GetDebug()>2) cerr<<"AliMUONDigitizer::Exec End of hits, track and file loops"<<endl;
333
334     // Loop on cathodes
335     Int_t icat;
336     for(icat=0; icat<2; icat++) {
337       //
338       // Filling Digit List
339       Int_t tracks[kMAXTRACKS];
340       Int_t charges[kMAXTRACKS];
341       Int_t nentries = fTDList->GetEntriesFast();
342       Int_t digits[6];
343       for (Int_t nent = 0; nent < nentries; nent++) {
344         AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
345         if (address == 0) continue; 
346         Int_t ich = address->Chamber();
347         Int_t   q = address->Signal(); 
348         chamber = &(pMUON->Chamber(ich));
349         //
350         //  Digit Response (noise, threshold, saturation, ...)
351         AliMUONResponse * response = chamber->ResponseModel();
352         q = response->DigitResponse(q,address);
353         
354         if (!q) continue;
355         
356         digits[0] = address->PadX();
357         digits[1] = address->PadY();
358         digits[2] = address->Cathode()-1;
359         digits[3] = q;
360         digits[4] = address->Physics();
361         digits[5] = address->Hit();
362         
363         Int_t nptracks = address->GetNTracks();
364         
365         if (nptracks > kMAXTRACKS) {
366           if (GetDebug() >0) {
367             cerr<<"AliMUONDigitizer:Exec  nptracks > 10 "<<nptracks;
368             cerr<<"reset to max value "<<kMAXTRACKS<<endl;
369           }
370           nptracks = kMAXTRACKS;
371         }
372         if (nptracks > 2 && GetDebug() >2) {
373           cerr<<"AliMUONDigitizer::Exec  nptracks > 2 "<<nptracks<<endl;
374           //    printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
375         }
376         for (Int_t tr = 0; tr < nptracks; tr++) {
377           tracks[tr]   = address->GetTrack(tr);
378           charges[tr]  = address->GetCharge(tr);
379         }      //end loop over list of tracks for one pad
380         // Sort list of tracks according to charge
381         if (nptracks > 1) {
382           SortTracks(tracks,charges,nptracks);
383         }
384         if (nptracks < kMAXTRACKS ) {
385           for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
386             tracks[i]  = -1;
387             charges[i] = 0;
388           }
389         }
390         
391         // fill digits
392         if (GetDebug()>3) cerr<<"AliMUONDigitzerv1::Exex TransientDigit to Digit"<<endl;
393         if ( digits[2] == icat ) pMUON->AddDigits(ich,tracks,charges,digits);
394       }
395       // Filling list of digits per chamber for a given cathode.
396       treeD->Fill();
397       pMUON->ResetDigits();    
398     } // end loop cathode
399     fTDList->Delete();  
400     
401     for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
402       if (fHitMap[ii]) {
403         delete fHitMap[ii];
404         fHitMap[ii] = 0;
405       }
406     }
407     
408     if (GetDebug()>2) 
409       cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
410           <<treeD->GetName()<<endl;
411
412     gimeout->WriteDigits("OVERWRITE");
413     delete [] fHitMap;
414     delete fTDList;
415     
416     if (fHits)    fHits->Clear();
417 }
418
419 //------------------------------------------------------------------------
420 void AliMUONDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
421 {
422   //
423   // Sort the list of tracks contributing to a given digit
424   // Only the 3 most significant tracks are acctually sorted
425   //
426   
427   //
428   //  Loop over signals, only 3 times
429   //
430   
431   Int_t qmax;
432   Int_t jmax;
433   Int_t idx[3] = {-2,-2,-2};
434   Int_t jch[3] = {-2,-2,-2};
435   Int_t jtr[3] = {-2,-2,-2};
436   Int_t i,j,imax;
437   
438   if (ntr<3) imax=ntr;
439   else imax=3;
440   for(i=0;i<imax;i++){
441     qmax=0;
442     jmax=0;
443     
444     for(j=0;j<ntr;j++){
445       
446       if((i == 1 && j == idx[i-1]) 
447          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
448       
449       if(charges[j] > qmax) {
450         qmax = charges[j];
451         jmax=j;
452       }       
453     } 
454     
455     if(qmax > 0) {
456       idx[i]=jmax;
457       jch[i]=charges[jmax]; 
458       jtr[i]=tracks[jmax];
459     }
460     
461   } 
462   
463   for(i=0;i<3;i++){
464     if (jtr[i] == -2) {
465       charges[i]=0;
466       tracks[i]=0;
467     } else {
468       charges[i]=jch[i];
469       tracks[i]=jtr[i];
470     }
471   }
472 }