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