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