Connect only MUON branches when reading the event to speed up digitisation.
[u/mrichter/AliRoot.git] / MUON / AliMUONDigitizer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.8  2002/02/13 09:03:24  jchudoba
19 Remove some deletes from dtor, those objects are deleted earlier in Exec() method (where they are created)
20
21 Revision 1.7  2001/11/22 11:15:41  jchudoba
22 Proper deletion of arrays (thanks to Rene Brun)
23
24 Revision 1.6  2001/11/02 12:55:45  jchudoba
25 cleanup of the code, add const to Get methods
26
27 Revision 1.4  2001/10/18 14:44:09  jchudoba
28 Define constant MAXTRACKS for maximum number of tracks associated with 1 digit
29
30 Revision 1.3  2001/10/04 20:01:54  jchudoba
31 changes for TTask implementation, some other small editing
32
33 Revision 1.2  2001/07/28 10:46:04  hristov
34 AliRunDigitizer.h included; typos corrected
35
36 Revision 1.1  2001/07/27 15:41:01  jchudoba
37 merging/digitization classes
38
39 */
40 #include <TTree.h> 
41 #include <TObjArray.h>
42 #include <TFile.h>
43 #include <TDirectory.h>
44 #include <iostream.h>
45
46 #include "AliMUONDigitizer.h"
47 #include "AliMUONConstants.h"
48 #include "AliMUONChamber.h"
49 #include "AliMUONHitMapA1.h"
50 #include "AliMUON.h"
51 #include "AliMUONHit.h"
52 #include "AliMUONPadHit.h"
53 #include "AliMUONDigit.h"
54 #include "AliMUONTransientDigit.h"
55 #include "AliRun.h"
56 #include "AliPDG.h"
57 #include "AliRunDigitizer.h"
58
59 ClassImp(AliMUONDigitizer)
60
61 //___________________________________________
62 AliMUONDigitizer::AliMUONDigitizer() :AliDigitizer()
63 {
64 // Default ctor - don't use it
65   fHits = 0;
66   fPadHits = 0;
67   fHitMap = 0;
68   fTDList = 0;
69 }
70
71 //___________________________________________
72 AliMUONDigitizer::AliMUONDigitizer(AliRunDigitizer* manager) 
73     :AliDigitizer(manager)
74 {
75 // ctor which should be used
76   fHitMap  = 0;
77   fTDList  = 0;
78   fHits    = 0;
79   fPadHits = 0;
80   fDebug   = 0; 
81   if (GetDebug()>2) 
82     cerr<<"AliMUONDigitizer::AliMUONDigitizer"
83         <<"(AliRunDigitizer* manager) was processed"<<endl;
84 }
85
86 //------------------------------------------------------------------------
87 AliMUONDigitizer::~AliMUONDigitizer()
88 {
89 // Destructor
90   if (fHits)       delete fHits;
91   if (fPadHits)    delete fPadHits;
92 }
93
94 //------------------------------------------------------------------------
95 Bool_t AliMUONDigitizer::Exists(const AliMUONPadHit *padhit) const
96 {
97   return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
98 }
99
100 //------------------------------------------------------------------------
101 void AliMUONDigitizer::Update(AliMUONPadHit *padhit)
102 {
103     AliMUONTransientDigit *pdigit = 
104       static_cast<AliMUONTransientDigit*>(
105       fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
106
107     // update charge
108     //
109     Int_t iqpad    = padhit->QPad();        // charge per pad
110     pdigit->AddSignal(iqpad);
111     pdigit->AddPhysicsSignal(iqpad);            
112
113     // update list of tracks
114     //
115     Int_t track, charge;    
116     track = fTrack+fMask;
117     if (fSignal) {
118       charge = iqpad;
119     } else {
120       charge = kBgTag;
121     }
122     pdigit->UpdateTrackList(track,charge);
123 }
124
125 //------------------------------------------------------------------------
126 void AliMUONDigitizer::CreateNew(AliMUONPadHit *padhit)
127 {
128 // Create new AliMUONTransientDigit and add it to the fTDList
129
130   fTDList->AddAtAndExpand(
131     new AliMUONTransientDigit(fNch,fDigits),fCounter);
132   fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
133   AliMUONTransientDigit* pdigit = 
134     (AliMUONTransientDigit*)fTDList->Last();
135   // list of tracks
136   Int_t track, charge;    
137   track = fTrack+fMask;
138   if (fSignal) {
139     charge = padhit->QPad();
140   } else {
141     charge = kBgTag;
142   }
143   pdigit->AddToTrackList(track,charge);
144   fCounter++;
145 }
146
147
148 //------------------------------------------------------------------------
149 Bool_t AliMUONDigitizer::Init()
150 {
151 // Initialization 
152     
153   fHits     = new TClonesArray("AliMUONHit",1000);
154   fPadHits  = new TClonesArray("AliMUONPadHit",1000);
155   return kTRUE;
156 }
157
158
159 //------------------------------------------------------------------------
160 //void AliMUONDigitizer::Digitize()
161 void AliMUONDigitizer::Exec(Option_t* option)
162 {
163
164   TString optionString = option;
165   if (optionString.Data() == "deb") {
166     cout<<"AliMUONDigitizer::Exec: called with option deb "<<endl;
167     fDebug = 3;
168   }
169   AliMUONChamber*   iChamber;
170   AliSegmentation*  segmentation;
171   
172   if (GetDebug()>2) cerr<<"   AliMUONDigitizer::Digitize() starts"<<endl;
173   fTDList = new TObjArray;
174
175   AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
176   if (!pMUON) {
177     cerr<<"AliMUONDigitizer::Digitize Error:"
178         <<" module MUON not found in the input file"<<endl;
179     return;
180   }
181   pMUON->MakeBranchInTreeD(fManager->GetTreeD());
182   fHitMap= new AliMUONHitMapA1* [AliMUONConstants::NCh()];
183
184   //
185   // loop over cathodes
186   //
187
188   for (int icat = 0; icat < 2; icat++) { 
189     fCounter = 0;
190     for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
191       iChamber = &(pMUON->Chamber(i));
192 //      if (!(iChamber->Nsec() == 1 && icat == 1)) {
193         segmentation = iChamber->SegmentationModel(icat+1);
194         fHitMap[i] = new AliMUONHitMapA1(segmentation, fTDList);
195 //      }
196     }
197
198
199 // Loop over files to digitize
200     fSignal = kTRUE;
201     for (Int_t inputFile=0; inputFile<fManager->GetNinputs();
202          inputFile++) {
203 // Connect MUON branches
204
205       if (inputFile > 0 ) fSignal = kFALSE;
206       TBranch *branchHits = 0;
207       TBranch *branchPadHits = 0;
208       TTree *treeH = fManager->GetInputTreeH(inputFile);
209       if (GetDebug()>2) {
210         cerr<<"   inputFile , cathode = "<<inputFile<<" "
211             <<icat<<endl;
212         cerr<<"   treeH, fHits "<<treeH<<" "<<fHits<<endl;
213       }
214       if (treeH && fHits) {
215         branchHits = treeH->GetBranch("MUON");
216         if (branchHits) {
217           fHits->Clear();
218           branchHits->SetAddress(&fHits);
219         }
220         else
221           Error("Exec","branch MUON was not found");
222       }
223       if (GetDebug()>2) cerr<<"   branchHits = "<<branchHits<<endl;
224
225       if (treeH && fPadHits) {
226         branchPadHits = treeH->GetBranch("MUONCluster");
227         if (branchPadHits) 
228           branchPadHits->SetAddress(&fPadHits);
229         else
230           Error("Exec","branch MUONCluster was not found");
231       }
232       if (GetDebug()>2) cerr<<"   branchPadHits = "<<branchPadHits<<endl;
233
234 //
235 //   Loop over tracks
236 //
237
238       Int_t ntracks = (Int_t) treeH->GetEntries();
239
240       for (fTrack = 0; fTrack < ntracks; fTrack++) {
241         if (GetDebug()>2) cerr<<"   fTrack = "<<fTrack<<endl;
242         fHits->Clear();
243         fPadHits->Clear();
244         branchHits->GetEntry(fTrack);
245         branchPadHits->GetEntry(fTrack);
246         
247 //
248 //   Loop over hits
249
250         AliMUONHit* mHit;
251         for(Int_t i = 0; i < fHits->GetEntriesFast(); ++i) {
252           mHit = static_cast<AliMUONHit*>(fHits->At(i));
253           fNch = mHit->Chamber()-1;  // chamber number
254           if (fNch > AliMUONConstants::NCh()-1) {
255             cerr<<"AliMUONDigitizer: ERROR: "
256                 <<"fNch > AliMUONConstants::NCh()-1, fNch, NCh(): "
257                 <<fNch<<", "<< AliMUONConstants::NCh()<<endl;
258             return;
259           }
260           iChamber = &(pMUON->Chamber(fNch));
261 //
262 // Loop over pad hits
263           for (AliMUONPadHit* mPad =
264                  (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
265                mPad;
266                mPad = (AliMUONPadHit*)pMUON->NextPad(fPadHits))
267           {
268             Int_t cathode  = mPad->Cathode();      // cathode number
269             Int_t ipx      = mPad->PadX();         // pad number on X
270             Int_t ipy      = mPad->PadY();         // pad number on Y
271             Int_t iqpad    = Int_t(mPad->QPad());  // charge per pad
272             if (cathode != (icat+1)) continue;
273
274             fMask = fManager->GetMask(inputFile);
275             fDigits[0] = ipx;
276             fDigits[1] = ipy;
277             if (!(fHitMap[fNch]->ValidateHit(fDigits[0], fDigits[1]))) continue;
278             fDigits[2] = icat;
279             fDigits[3] = iqpad;
280             if (inputFile == 0) {
281               fDigits[4] = iqpad;
282             } else {
283               fDigits[4] = 0;
284             }
285             if (mHit->Particle() == kMuonPlus ||
286                 mHit->Particle() == kMuonMinus) {
287               fDigits[5] = (mPad->HitNumber()) + fMask; 
288             } else fDigits[5] = -1;
289
290             // build the list of fired pads and update the info, 
291             // fDigits is input for Update(mPad)
292
293             if (!Exists(mPad)) {
294               CreateNew(mPad);
295             } else {
296               Update(mPad);
297             } //  end if Exists(mPad)
298           } //end loop over clusters
299         } // hit loop
300       } // track loop
301     } // end file loop
302     if (GetDebug()>2) cerr<<"END OF FILE LOOP"<<endl;
303
304     Int_t tracks[kMAXTRACKS];
305     Int_t charges[kMAXTRACKS];
306     Int_t nentries = fTDList->GetEntriesFast();
307         
308     for (Int_t nent = 0; nent < nentries; nent++) {
309       AliMUONTransientDigit *address = (AliMUONTransientDigit*)fTDList->At(nent);
310       if (address == 0) continue; 
311       Int_t ich = address->Chamber();
312       Int_t   q = address->Signal(); 
313       iChamber = &(pMUON->Chamber(ich));
314 //
315 //  Digit Response (noise, threshold, saturation, ...)
316       AliMUONResponse * response = iChamber->ResponseModel();
317       q = response->DigitResponse(q);
318             
319       if (!q) continue;
320             
321       fDigits[0] = address->PadX();
322       fDigits[1] = address->PadY();
323       fDigits[2] = address->Cathode();
324       fDigits[3] = q;
325       fDigits[4] = address->Physics();
326       fDigits[5] = address->Hit();
327             
328       Int_t nptracks = address->GetNTracks();
329
330       if (nptracks > kMAXTRACKS) {
331         if (GetDebug() >0) {
332           cerr<<"AliMUONDigitizer: nptracks > 10 "<<nptracks;
333           cerr<<"reset to max value "<<kMAXTRACKS<<endl;
334         }
335         nptracks = kMAXTRACKS;
336       }
337       if (nptracks > 2 && GetDebug() >2) {
338         cerr<<"AliMUONDigitizer: nptracks > 2 "<<nptracks<<endl;
339         printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,fDigits[0],fDigits[1],q);
340       }
341       for (Int_t tr = 0; tr < nptracks; tr++) {
342         tracks[tr]   = address->GetTrack(tr);
343         charges[tr]  = address->GetCharge(tr);
344       }      //end loop over list of tracks for one pad
345       // Sort list of tracks according to charge
346       if (nptracks > 1) {
347         SortTracks(tracks,charges,nptracks);
348       }
349       if (nptracks < kMAXTRACKS ) {
350         for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
351           tracks[i]  = 0;
352           charges[i] = 0;
353         }
354       }
355             
356       // fill digits
357       pMUON->AddDigits(ich,tracks,charges,fDigits);
358     }
359
360     fManager->GetTreeD()->Fill();
361
362     pMUON->ResetDigits();  //
363     fTDList->Clear();
364
365         
366     for(Int_t ii = 0; ii < AliMUONConstants::NCh(); ++ii) {
367       if (fHitMap[ii]) {
368         delete fHitMap[ii];
369         fHitMap[ii] = 0;
370       }
371     }
372   } //end loop over cathodes
373   if (GetDebug()>2) 
374     cerr<<"AliMUONDigitizer::Exec: writing the TreeD: "
375         <<fManager->GetTreeD()->GetName()<<endl;
376   fManager->GetTreeD()->Write(0,TObject::kOverwrite);
377   delete [] fHitMap;
378   delete fTDList;
379     
380   if (fHits)    fHits->Delete();
381   if (fPadHits) fPadHits->Delete();
382 }
383
384
385 //------------------------------------------------------------------------
386 void AliMUONDigitizer::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
387 {
388   //
389   // Sort the list of tracks contributing to a given digit
390   // Only the 3 most significant tracks are acctually sorted
391   //
392   
393   //
394   //  Loop over signals, only 3 times
395   //
396   
397   Int_t qmax;
398   Int_t jmax;
399   Int_t idx[3] = {-2,-2,-2};
400   Int_t jch[3] = {-2,-2,-2};
401   Int_t jtr[3] = {-2,-2,-2};
402   Int_t i,j,imax;
403   
404   if (ntr<3) imax=ntr;
405   else imax=3;
406   for(i=0;i<imax;i++){
407     qmax=0;
408     jmax=0;
409     
410     for(j=0;j<ntr;j++){
411       
412       if((i == 1 && j == idx[i-1]) 
413          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
414       
415       if(charges[j] > qmax) {
416         qmax = charges[j];
417         jmax=j;
418       }       
419     } 
420     
421     if(qmax > 0) {
422       idx[i]=jmax;
423       jch[i]=charges[jmax]; 
424       jtr[i]=tracks[jmax]; 
425     }
426     
427   } 
428   
429   for(i=0;i<3;i++){
430     if (jtr[i] == -2) {
431       charges[i]=0;
432       tracks[i]=0;
433     } else {
434       charges[i]=jch[i];
435       tracks[i]=jtr[i];
436     }
437   }
438 }