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