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