Bug on digit output. Now treeD event per cathode
[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->Clear();
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     // Loop on cathodes
276     Int_t icat;
277     for(icat=0; icat<2; icat++) {
278       //
279       // Filling Digit List
280       Int_t tracks[kMAXTRACKS];
281       Int_t charges[kMAXTRACKS];
282       Int_t nentries = fTDList->GetEntriesFast();
283       Int_t digits[6];
284       for (Int_t nent = 0; nent < nentries; nent++) {
285         AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
286         if (address == 0) continue; 
287         Int_t ich = address->Chamber();
288         Int_t   q = address->Signal(); 
289         chamber = &(pMUON->Chamber(ich));
290         //
291         //  Digit Response (noise, threshold, saturation, ...)
292         AliMUONResponse * response = chamber->ResponseModel();
293         q = response->DigitResponse(q,address);
294         
295         if (!q) continue;
296         
297         digits[0] = address->PadX();
298         digits[1] = address->PadY();
299         digits[2] = address->Cathode()-1;
300         digits[3] = q;
301         digits[4] = address->Physics();
302         digits[5] = address->Hit();
303         
304         Int_t nptracks = address->GetNTracks();
305         
306         if (nptracks > kMAXTRACKS) {
307           if (GetDebug() >0) {
308             cerr<<"AliMUONDigitizer:Exec  nptracks > 10 "<<nptracks;
309             cerr<<"reset to max value "<<kMAXTRACKS<<endl;
310           }
311           nptracks = kMAXTRACKS;
312         }
313         if (nptracks > 2 && GetDebug() >2) {
314           cerr<<"AliMUONDigitizer::Exec  nptracks > 2 "<<nptracks<<endl;
315           //    printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
316         }
317         for (Int_t tr = 0; tr < nptracks; tr++) {
318           tracks[tr]   = address->GetTrack(tr);
319           charges[tr]  = address->GetCharge(tr);
320         }      //end loop over list of tracks for one pad
321         // Sort list of tracks according to charge
322         if (nptracks > 1) {
323           SortTracks(tracks,charges,nptracks);
324         }
325         if (nptracks < kMAXTRACKS ) {
326           for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
327             tracks[i]  = 0;
328             charges[i] = 0;
329           }
330         }
331         
332         // fill digits
333         if (GetDebug()>2) cerr<<"AliMUONDigitzerv1::Exex TransientDigit to Digit"<<endl;
334         if ( digits[2] == icat ) pMUON->AddDigits(ich,tracks,charges,digits);
335       }
336       fManager->GetTreeD()->Fill();
337       pMUON->ResetDigits();  //   
338     } // end loop cathode
339     fTDList->Delete();  
340     
341     for(Int_t ii = 0; ii < 2*AliMUONConstants::NCh(); ++ii) {
342       if (fHitMap[ii]) {
343         delete fHitMap[ii];
344         fHitMap[ii] = 0;
345       }
346     }
347     
348     if (GetDebug()>2) 
349       cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
350           <<fManager->GetTreeD()->GetName()<<endl;
351     fManager->GetTreeD()->Write(0,TObject::kOverwrite);
352     delete [] fHitMap;
353     delete fTDList;
354     
355     if (fHits)    fHits->Clear();
356 }
357
358 //------------------------------------------------------------------------
359 void AliMUONDigitizerv1::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
360 {
361   //
362   // Sort the list of tracks contributing to a given digit
363   // Only the 3 most significant tracks are acctually sorted
364   //
365   
366   //
367   //  Loop over signals, only 3 times
368   //
369   
370   Int_t qmax;
371   Int_t jmax;
372   Int_t idx[3] = {-2,-2,-2};
373   Int_t jch[3] = {-2,-2,-2};
374   Int_t jtr[3] = {-2,-2,-2};
375   Int_t i,j,imax;
376   
377   if (ntr<3) imax=ntr;
378   else imax=3;
379   for(i=0;i<imax;i++){
380     qmax=0;
381     jmax=0;
382     
383     for(j=0;j<ntr;j++){
384       
385       if((i == 1 && j == idx[i-1]) 
386          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
387       
388       if(charges[j] > qmax) {
389         qmax = charges[j];
390         jmax=j;
391       }       
392     } 
393     
394     if(qmax > 0) {
395       idx[i]=jmax;
396       jch[i]=charges[jmax]; 
397       jtr[i]=tracks[jmax]; 
398     }
399     
400   } 
401   
402   for(i=0;i<3;i++){
403     if (jtr[i] == -2) {
404       charges[i]=0;
405       tracks[i]=0;
406     } else {
407       charges[i]=jch[i];
408       tracks[i]=jtr[i];
409     }
410   }
411 }