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