e0dd5a16328a87d8c634b7858ed3b2c7b2ad7b85
[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
10 #include "AliMUON.h"
11 #include "AliMUONChamber.h"
12 #include "AliMUONConstants.h"
13 #include "AliMUONDigit.h"
14 #include "AliMUONDigitizerv1.h"
15 #include "AliMUONHit.h"
16 #include "AliMUONHitMapA1.h"
17 #include "AliMUONPadHit.h"
18 #include "AliMUONTransientDigit.h"
19 #include "AliRun.h"
20 #include "AliRunDigitizer.h"
21
22 ClassImp(AliMUONDigitizerv1)
23
24 //___________________________________________
25 AliMUONDigitizerv1::AliMUONDigitizerv1() :AliDigitizer()
26 {
27 // Default ctor - don't use it
28   fHitMap = 0;
29   fTDList = 0;
30   if (GetDebug()>2) 
31     cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
32         <<"(AliRunDigitizer* manager) was processed"<<endl;
33 }
34
35 //___________________________________________
36 AliMUONDigitizerv1::AliMUONDigitizerv1(AliRunDigitizer* manager) 
37     :AliDigitizer(manager)
38 {
39 // ctor which should be used
40   fHitMap  = 0;
41   fTDList  = 0;
42   fDebug   = 0; 
43   fHits = new TClonesArray("AliMUONHit",1000);
44   if (GetDebug()>2) 
45     cerr<<"AliMUONDigitizerv1::AliMUONDigitizerv1"
46         <<"(AliRunDigitizer* manager) was processed"<<endl;
47 }
48
49 //------------------------------------------------------------------------
50 AliMUONDigitizerv1::~AliMUONDigitizerv1()
51 {
52 // Destructor
53   delete fHits;
54 }
55
56 //------------------------------------------------------------------------
57 void AliMUONDigitizerv1::AddTransientDigit(AliMUONTransientDigit * mTD)
58 {
59   // Choosing the maping of the cathode plane of the chamber:
60   Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
61   fTDList->AddAtAndExpand(mTD, fTDCounter);
62   fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
63   fTDCounter++;
64 }
65
66 //------------------------------------------------------------------------
67 Bool_t AliMUONDigitizerv1::ExistTransientDigit(AliMUONTransientDigit * mTD)
68 {
69   // Choosing the maping of the cathode plane of the chamber:
70   Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
71   return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
72 }
73
74 //------------------------------------------------------------------------
75 Bool_t AliMUONDigitizerv1::Init()
76 {
77 // Initialization 
78   return kTRUE;
79 }
80
81 //------------------------------------------------------------------------
82 void AliMUONDigitizerv1::UpdateTransientDigit(Int_t track, AliMUONTransientDigit * mTD)
83 {
84   // Choosing the maping of the cathode plane of the chamber:
85   Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
86   AliMUONTransientDigit *pdigit = 
87     static_cast<AliMUONTransientDigit*>(fHitMap[iNchCpl]->GetHit(mTD->PadX(),mTD->PadY()));
88   // update charge
89   //
90   Int_t iqpad  = mTD->Signal();        // charge per pad
91   pdigit->AddSignal(iqpad);
92   pdigit->AddPhysicsSignal(iqpad);              
93   // update list of tracks
94   //
95   Int_t charge;    
96   track=+ fMask;
97   if (fSignal) charge = iqpad;
98   else         charge = kBgTag;
99   pdigit->UpdateTrackList(track,charge);
100 }
101
102
103 //--------------------------------------------------------------------------
104 void AliMUONDigitizerv1::MakeTransientDigit(Int_t track, Int_t iHit, AliMUONHit * mHit)
105 {
106   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
107   if (!pMUON) { 
108     cerr<<"AliMUONDigitizerv1::Digitize Error:"
109         <<" module MUON not found in the input file"<<endl;
110   }
111   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit starts"<<endl;
112   Int_t   ichamber = mHit->Chamber()-1;
113   AliMUONChamber & chamber = pMUON->Chamber(ichamber);
114   Float_t xhit = mHit->X();
115   Float_t yhit = mHit->Y();
116   Float_t zhit = mHit->Z();
117   Float_t eloss= mHit->Eloss(); 
118   Float_t tof  = mHit->Age();
119   // Variables for chamber response from AliMUONChamber::DisIntegration
120   Float_t newdigit[6][500];  // Pad information
121   Int_t nnew=0;              // Number of touched Pads per hit
122   Int_t digits[6];
123   
124   //
125   // Calls the charge disintegration method of the current chamber 
126   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit calling AliMUONChamber::DisIngtegration starts"<<endl;
127   chamber.DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newdigit);
128   // Creating a new TransientDigits from hit
129   for(Int_t iTD=0; iTD<nnew; iTD++) {
130     digits[0] = Int_t(newdigit[1][iTD]);  // Padx of the Digit
131     digits[1] = Int_t(newdigit[2][iTD]);  // Pady of the Digit
132     digits[2] = Int_t(newdigit[5][iTD]);  // Cathode plane
133     digits[3] = Int_t(newdigit[3][iTD]);  // Induced charge in the Pad
134     if (fSignal) digits[4] = Int_t(newdigit[3][iTD]);
135     else         digits[4] = 0;
136     digits[5] = iHit+fMask;    // Hit number in the list
137     if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit " <<
138                         "PadX "<< digits[0] << " " <<
139                         "PadY "<< digits[1] << " " <<
140                         "Plane " << digits[2] << " " <<
141                         "Charge " << digits[3] <<" " <<
142                         "Hit " << digits[5] << endl;
143     // list of tracks
144     Int_t charge;    
145     track += fMask;
146     if (fSignal) charge = digits[3];
147     else         charge = kBgTag;
148     if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Creating AliMUONTransientDigit"<<endl;
149     AliMUONTransientDigit * mTD = new AliMUONTransientDigit(ichamber, digits);
150     mTD->AddToTrackList(track,charge);
151
152     if (!ExistTransientDigit(mTD)) {
153       AddTransientDigit(mTD);
154       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit Adding TransientDigit"<<endl;
155     }
156     else {
157       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::MakeTransientDigit updating TransientDigit"<<endl;
158       UpdateTransientDigit(track, mTD);
159       delete mTD;
160     }
161   }
162 }
163 //-----------------------------------------------------------------------
164 void AliMUONDigitizerv1::Exec(Option_t* option)
165 {
166   TString optionString = option;
167   if (optionString.Data() == "deb") {
168     cout<<"AliMUONDigitizerv1::Exec: called with option deb "<<endl;
169     fDebug = 3;
170   }
171
172   AliMUONChamber*   chamber;
173   AliSegmentation*  c1Segmentation; //Cathode plane c1 of the chamber
174   AliSegmentation*  c2Segmentation; //Cathode place c2 of the chamber
175
176   if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Digitize() starts"<<endl;
177   fTDList = new TObjArray;
178
179   // Getting Module MUON
180   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
181   if (!pMUON) {
182     cerr<<"AliMUONDigitizerv1::Digitize Error:"
183         <<" module MUON not found in the input file"<<endl;
184     return;
185   }
186   // New branch for MUON digit in the tree of digits
187   pMUON->MakeBranchInTreeD(fManager->GetTreeD());
188
189   // Array of pointer of the AliMUONHitMapA1:
190   //  two HitMaps per chamber, or one HitMap per cahtode plane
191   fHitMap= new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];
192
193   //Loop over chambers for the definition AliMUONHitMap
194   for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
195     chamber = &(pMUON->Chamber(i));
196     c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
197     fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
198     c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
199     fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
200   }
201
202 // Loop over files to merge and to digitize
203     fSignal = kTRUE;
204     for (Int_t inputFile=0; inputFile<fManager->GetNinputs(); inputFile++) {
205       // Connect MUON Hit branch
206       if (inputFile > 0 ) fSignal = kFALSE;
207       TBranch *branchHits = 0;
208       TTree *treeH = fManager->GetInputTreeH(inputFile);
209       if (GetDebug()>2) {
210         cerr<<"AliMUONDigitizerv1::Exec inputFile is "<<inputFile<<" "<<endl;
211         cerr<<"AliMUONDigitizerv1::Exec treeH, fHits "<<treeH<<" "<<fHits<<endl;
212       }
213       if (treeH && fHits) {
214         branchHits = treeH->GetBranch("MUON");
215         if (branchHits) {
216           fHits->Delete();
217           branchHits->SetAddress(&fHits);
218         }
219         else
220           Error("Exec","branch MUON was not found");
221       }
222       if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec branchHits = "<<branchHits<<endl;
223
224       fMask = fManager->GetMask(inputFile);
225       //
226       // Loop over tracks
227       Int_t itrack;
228       Int_t ntracks = (Int_t) treeH->GetEntries();
229       for (itrack = 0; itrack < ntracks; itrack++) {
230         if (GetDebug()>2) cerr<<"AliMUONDigitizerv1::Exec itrack = "<<itrack<<endl;
231         fHits->Clear();
232         branchHits->GetEntry(itrack);
233         //
234         //  Loop over hits
235         Int_t ihit, ichamber;
236         AliMUONHit* mHit;
237         for(ihit = 0; ihit < fHits->GetEntriesFast(); ihit++) {
238           mHit = static_cast<AliMUONHit*>(fHits->At(ihit));
239           ichamber = mHit->Chamber()-1;  // chamber number
240           if (ichamber > AliMUONConstants::NCh()-1) {
241             cerr<<"AliMUONDigitizer: ERROR: "
242                 <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
243                 <<ichamber<<", "<< AliMUONConstants::NCh()<<endl;
244             return;
245           }
246           chamber = &(pMUON->Chamber(ichamber));
247           //
248           //Dumping Hit content:
249           if (GetDebug()>2) {
250             cerr<<"AliMuonDigitizerv1::Exec ihit, ichamber, x, y, z, eloss " <<
251               ihit << " " << 
252               mHit->Chamber() << " " <<
253               mHit->X() << " " <<
254               mHit->Y() << " " <<
255               mHit->Z() << " " <<
256               mHit->Eloss() << " " << endl;
257           }
258           // 
259           // Inititializing Correlation
260           chamber->ChargeCorrelationInit();
261           if (ichamber < AliMUONConstants::NTrackingCh()) {
262             // Tracking Chamber
263             // Initialize hit position (cursor) in the segmentation model 
264             chamber->SigGenInit(mHit->X(), mHit->Y(), mHit->Z());
265           } else {
266             // Trigger Chamber 
267           }
268           MakeTransientDigit(itrack, ihit, mHit);
269         } // hit loop
270       } // track loop
271     } // end file loop
272     if (GetDebug()>2) cerr<<"AliMUONDigitizer::Exec End of hits, track and file loops"<<endl;
273
274     //
275     // Filling Digit List
276     Int_t tracks[kMAXTRACKS];
277     Int_t charges[kMAXTRACKS];
278     Int_t nentries = fTDList->GetEntriesFast();
279     Int_t digits[6];
280     for (Int_t nent = 0; nent < nentries; nent++) {
281       AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
282       if (address == 0) continue; 
283       Int_t ich = address->Chamber();
284       Int_t   q = address->Signal(); 
285       chamber = &(pMUON->Chamber(ich));
286       //
287       //  Digit Response (noise, threshold, saturation, ...)
288       AliMUONResponse * response = chamber->ResponseModel();
289       q = response->DigitResponse(q,address);
290       
291       if (!q) continue;
292       
293       digits[0] = address->PadX();
294       digits[1] = address->PadY();
295       digits[2] = address->Cathode();
296       digits[3] = q;
297       digits[4] = address->Physics();
298       digits[5] = address->Hit();
299       
300       Int_t nptracks = address->GetNTracks();
301       
302       if (nptracks > kMAXTRACKS) {
303         if (GetDebug() >0) {
304           cerr<<"AliMUONDigitizer:Exec  nptracks > 10 "<<nptracks;
305           cerr<<"reset to max value "<<kMAXTRACKS<<endl;
306         }
307         nptracks = kMAXTRACKS;
308       }
309       if (nptracks > 2 && GetDebug() >2) {
310         cerr<<"AliMUONDigitizer::Exec  nptracks > 2 "<<nptracks<<endl;
311         //      printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
312       }
313       for (Int_t tr = 0; tr < nptracks; tr++) {
314         tracks[tr]   = address->GetTrack(tr);
315         charges[tr]  = address->GetCharge(tr);
316       }      //end loop over list of tracks for one pad
317       // Sort list of tracks according to charge
318       if (nptracks > 1) {
319         SortTracks(tracks,charges,nptracks);
320       }
321       if (nptracks < kMAXTRACKS ) {
322         for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
323           tracks[i]  = 0;
324           charges[i] = 0;
325         }
326       }
327       
328       // fill digits
329       if (GetDebug()>2) cerr<<"AliMUONDigitzerv1::Exex TransientDigit to Digit"<<endl;
330       pMUON->AddDigits(ich,tracks,charges,digits);
331     }
332     fManager->GetTreeD()->Fill();
333     pMUON->ResetDigits();  //
334     fTDList->Delete();
335     
336     for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
337       if (fHitMap[ii]) {
338         delete fHitMap[ii];
339         fHitMap[ii] = 0;
340       }
341     }
342     
343     if (GetDebug()>2) 
344       cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
345           <<fManager->GetTreeD()->GetName()<<endl;
346     fManager->GetTreeD()->Write(0,TObject::kOverwrite);
347     delete [] fHitMap;
348     delete fTDList;
349     
350     if (fHits)    fHits->Clear();
351 }
352
353 //------------------------------------------------------------------------
354 void AliMUONDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
355 {
356   //
357   // Sort the list of tracks contributing to a given digit
358   // Only the 3 most significant tracks are acctually sorted
359   //
360   
361   //
362   //  Loop over signals, only 3 times
363   //
364   
365   Int_t qmax;
366   Int_t jmax;
367   Int_t idx[3] = {-2,-2,-2};
368   Int_t jch[3] = {-2,-2,-2};
369   Int_t jtr[3] = {-2,-2,-2};
370   Int_t i,j,imax;
371   
372   if (ntr<3) imax=ntr;
373   else imax=3;
374   for(i=0;i<imax;i++){
375     qmax=0;
376     jmax=0;
377     
378     for(j=0;j<ntr;j++){
379       
380       if((i == 1 && j == idx[i-1]) 
381          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
382       
383       if(charges[j] > qmax) {
384         qmax = charges[j];
385         jmax=j;
386       }       
387     } 
388     
389     if(qmax > 0) {
390       idx[i]=jmax;
391       jch[i]=charges[jmax]; 
392       jtr[i]=tracks[jmax]; 
393     }
394     
395   } 
396   
397   for(i=0;i<3;i++){
398     if (jtr[i] == -2) {
399       charges[i]=0;
400       tracks[i]=0;
401     } else {
402       charges[i]=jch[i];
403       tracks[i]=jtr[i];
404     }
405   }
406 }