Got rid of class template AliFMD<Type> on request of Federico, who
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 //____________________________________________________________________
19 //
20 // This is a class that constructs ReconstParticles (reconstructed
21 // particles) out of Digits
22 //
23 //
24 // This class reads either digits from a TClonesArray or raw data from
25 // a DDL file (or similar), and stores the read ADC counts in an
26 // internal cache (fAdcs). 
27 //
28 // From the cached values it then calculates the number of particles
29 // that hit a region of the FMDs, as specified by the user. 
30 //
31 // The reconstruction can be done in two ways: Either via counting the
32 // number of empty strips (Poisson method), or by converting the ADC
33 // signal to an energy deposition, and then dividing by the typical
34 // energy loss of a particle.
35 // 
36 // Currently, this class only reads the digits from a TClonesArray,
37 // and the Poission method for reconstruction. 
38 //
39 // 
40 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
41 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
42 //
43 //
44 //____________________________________________________________________
45
46 #include "AliFMD.h"                             // ALIFMD_H
47 #include "AliFMDDigit.h"                        // ALIFMDDIGIT_H
48 #include "AliFMDParticles.h"                    // ALIFMDPARTICLES_H
49 #include "AliFMDReconstructor.h"                // ALIFMDRECONSTRUCTOR_H
50 #include "AliAltroBuffer.h"                     // ALIALTROBUFFER_H
51 #include "AliLog.h"                             // ALILOG_H
52 #include "AliRun.h"                             // ALIRUN_H
53 #include "AliRunLoader.h"                       // ALIRUNLOADER_H
54 #include "AliLoader.h"                          // ALILOADER_H
55 #include "AliHeader.h"                          // ALIHEADER_H
56 #include "AliGenEventHeader.h"                  // ALIGENEVENTHEADER_H
57 #include "AliFMDRawStream.h"                    // ALIFMDRAWSTREAM_H
58 #include "AliFMDRawReader.h"                    // ALIFMDRAWREADER_H
59 #include "AliRawReader.h"                       // ALIRAWREADER_H
60 #include "AliFMDReconstructionAlgorithm.h"      // ALIFMDRECONSTRUCTIONALGORITHM_H
61
62 //____________________________________________________________________
63 ClassImp(AliFMDReconstructor);
64
65 //____________________________________________________________________
66 AliFMDReconstructor::AliFMDReconstructor() 
67   : AliReconstructor(),
68     fDeltaEta(0), 
69     fDeltaPhi(0), 
70     fThreshold(0),
71     fPedestal(0), 
72     fPedestalWidth(0),
73     fPedestalFactor(0)
74 {
75   // Make a new FMD reconstructor object - default CTOR.
76   SetDeltaEta();
77   SetDeltaPhi();
78   SetThreshold();
79   SetPedestal();
80
81   fParticles = new TClonesArray("AliFMDParticles", 1000);
82   fFMDLoader = 0;
83   fRunLoader = 0;
84   fFMD       = 0;
85 }
86   
87
88 //____________________________________________________________________
89 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
90   : AliReconstructor(),
91     fDeltaEta(0), 
92     fDeltaPhi(0), 
93     fThreshold(0),
94     fPedestal(0), 
95     fPedestalWidth(0),
96     fPedestalFactor(0)
97 {
98   // Make a new FMD reconstructor object - default CTOR.
99   SetDeltaEta(other.fDeltaEta);
100   SetDeltaPhi(other.fDeltaPhi);
101   SetThreshold(other.fThreshold);
102   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
103
104   // fParticles = new TClonesArray("AliFMDParticles", 1000);
105   fFMDLoader = other.fFMDLoader;
106   fRunLoader = other.fRunLoader;
107   fFMD       = other.fFMD;
108 }
109   
110
111 //____________________________________________________________________
112 AliFMDReconstructor&
113 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
114 {
115   // Make a new FMD reconstructor object - default CTOR.
116   SetDeltaEta(other.fDeltaEta);
117   SetDeltaPhi(other.fDeltaPhi);
118   SetThreshold(other.fThreshold);
119   SetPedestal(other.fPedestal, other.fPedestalWidth);
120
121   // fParticles = new TClonesArray("AliFMDParticles", 1000);
122   fFMDLoader = other.fFMDLoader;
123   fRunLoader = other.fRunLoader;
124   fFMD       = other.fFMD;
125
126   return *this;
127 }
128   
129 //____________________________________________________________________
130 void 
131 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor) 
132 {
133   // Set the pedestal, and pedestal width 
134   fPedestal       = mean;
135   fPedestalWidth  = width;
136   fPedestalFactor = factor;
137 }
138
139 //____________________________________________________________________
140 void 
141 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader, 
142                                  AliRawReader* rawReader) const
143
144   // Collects all digits in the same active volume into number of
145   // particles
146   //
147   // Reconstruct number of particles in given group of pads for given
148   // FMDvolume determined by numberOfVolume,
149   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
150   // numberOgMaxRing 
151   //
152   // The reconstruction method is choosen based on the number of empty
153   // strips. 
154   fParticles->Clear();
155   if (!runLoader) {
156     Error("Exec","Run Loader loader is NULL - Session not opened");
157     return;
158   }
159   fRunLoader = runLoader;
160   fFMDLoader = runLoader->GetLoader("FMDLoader");
161   if (!fFMDLoader) 
162     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
163           "in specified event");
164
165   // Get the AliRun object
166   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
167   
168   // Get the AliFMD object
169   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
170   if (!fFMD) {
171     AliError("Can not get FMD from gAlice");
172     return;
173   }
174   fFMDLoader->LoadRecPoints("RECREATE");
175
176   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
177
178   if (rawReader) {
179     Int_t event = 0;
180     while (rawReader->NextEvent()) {
181       ProcessEvent(event, rawReader);
182       event++;
183     }
184   }
185   else {
186     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries()); 
187     for(Int_t event = 0; event < nEvents; event++) 
188       ProcessEvent(event, 0);
189   }
190
191
192   fFMDLoader->UnloadRecPoints();
193   fFMDLoader = 0;
194   fRunLoader = 0;
195   fFMD       = 0;
196 }
197
198 //____________________________________________________________________
199 void 
200 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
201
202   // Collects all digits in the same active volume into number of
203   // particles
204   //
205   // Reconstruct number of particles in given group of pads for given
206   // FMDvolume determined by numberOfVolume,
207   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
208   // numberOgMaxRing 
209   //
210   // The reconstruction method is choosen based on the number of empty
211   // strips. 
212   Reconstruct(runLoader, 0);
213 }
214
215
216 //____________________________________________________________________
217 void 
218 AliFMDReconstructor::ProcessEvent(Int_t event, 
219                                   AliRawReader* reader) const
220 {
221   // Process one event read from either a clones array or from a a raw
222   // data reader. 
223   fRunLoader->GetEvent(event) ;
224   //event z-vertex for correction eta-rad dependence      
225   AliHeader *header            = fRunLoader->GetHeader();
226   if (!header) Warning("ProcessEvent", "no AliHeader found!");
227   AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
228
229   // Get the Z--coordinate from the event header 
230   TArrayF o(3); 
231   if (genHeader) genHeader->PrimaryVertex(o);
232   else           Warning("ProcessEvent", "No GenEventHeader Found");
233   fCurrentVertex = o.At(2);
234
235   // If the recontruction tree isn't loaded, load it
236   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
237   
238   //Make branches to hold the reconstructed particles 
239   const Int_t kBufferSize = 16000;
240   fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
241
242   // Load or recreate the digits 
243   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
244     if (!reader) {
245       Error("Exec","Error occured while loading digits. Exiting.");
246       return;
247     }
248   }
249   // Get the digits tree 
250   TTree* digitTree = fFMDLoader->TreeD();
251   if (!digitTree) { 
252     if (!reader) {
253       Error("Exec","Can not get Tree with Digits. "
254             "Nothing to reconstruct - Exiting");
255       return;
256     }
257     fFMDLoader->MakeTree("D");
258     digitTree = fFMDLoader->TreeD();
259     
260   }
261   // Get the FMD branch holding the digits. 
262   TBranch *digitBranch = digitTree->GetBranch("FMD");
263   TClonesArray* digits = fFMD->Digits();
264   if (!digitBranch) {
265     if (!reader) {
266       Error("Exec", "No digit branch for the FMD found");
267       return;
268     }
269     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
270   }
271   if (!reader) digitBranch->SetAddress(&digits);
272
273   if  (reader) {
274     AliFMDRawReader rawRead(fFMD, reader);
275     // rawRead->SetSampleRate(fSampleRate);
276     rawRead.Exec();
277   }
278   else {
279     // Read the ADC values from a clones array. 
280     AliDebug(10, "Reading ADCs from Digits array");
281     // read Digits, and reconstruct the particles
282     if (!fFMDLoader->TreeD()->GetEvent(0)) return;
283   }
284   
285   TIter next(&fAlgorithms);
286   AliFMDReconstructionAlgorithm* algorithm = 0;
287   while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next()))) 
288     algorithm->PreEvent();
289
290   ProcessDigits(digits);
291
292   next.Reset();
293   algorithm = 0;
294   while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next()))) 
295     algorithm->PostEvent();
296   
297   if (reader) {
298     digitTree->Fill();
299     fFMDLoader->WriteDigits("OVERWRITE");
300   }
301   fFMDLoader->UnloadDigits();
302   fFMDLoader->TreeR()->Reset();
303   fFMDLoader->TreeR()->Fill(); 
304   fFMDLoader->WriteRecPoints("OVERWRITE");
305 }
306
307 //____________________________________________________________________
308 UShort_t
309 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
310 {
311   // Member function to subtract the pedestal from a digit
312   // This implementation does nothing, but a derived class could over
313   // load this to subtract a pedestal that was given in a database or
314   // something like that. 
315
316   Int_t counts = 0;
317   Float_t ped = fPedestal * fPedestalFactor * fPedestalWidth;
318   if (digit->Count3() >= 0)      counts = digit->Count3();
319   else if (digit->Count2() >= 0) counts = digit->Count2();
320   else                           counts = digit->Count2();
321   counts = TMath::Max(Int_t(counts - ped), 0);
322   return  UShort_t(counts);
323 }
324
325 //____________________________________________________________________
326 void
327 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
328 {
329   Int_t nDigits = digits->GetEntries();
330   for (Int_t i = 0; i < nDigits; i++) {
331     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
332     AliFMDSubDetector* subDetector = 0;
333     switch (digit->Detector()) {
334     case 1: subDetector = fFMD->GetFMD1(); break;
335     case 2: subDetector = fFMD->GetFMD2(); break;
336     case 3: subDetector = fFMD->GetFMD3(); break;
337     }
338     if (!subDetector) { 
339       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
340       continue;
341     }
342     
343     AliFMDRing* ring  = 0;
344     Float_t     ringZ = 0;
345     switch(digit->Ring()) {
346     case 'i':
347     case 'I':  
348       ring  = subDetector->GetInner(); 
349       ringZ = subDetector->GetInnerZ();
350       break;
351     case 'o':
352     case 'O':  
353       ring  = subDetector->GetOuter(); 
354       ringZ = subDetector->GetOuterZ();
355       break;
356     }
357     if (!ring) {
358       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
359               digit->Ring());
360       break;
361     }
362     
363     Float_t  realZ    = fCurrentVertex + ringZ;
364     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) / ring->GetNStrips() 
365                          * (digit->Strip() + .5) + ring->GetLowR());
366     Float_t  theta    = TMath::ATan2(stripR, realZ);
367     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
368                          * (digit->Sector() + .5));
369     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
370     UShort_t counts   = SubtractPedestal(digit);
371     
372     TIter next(&fAlgorithms);
373     AliFMDReconstructionAlgorithm* algorithm = 0;
374     while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next()))) 
375       algorithm->ProcessDigit(digit, eta, phi, counts);
376   }
377 }
378
379       
380 //____________________________________________________________________
381 void
382 AliFMDReconstructor::ReconstructFromCache() const
383 {
384   // Based on the information in the cache, do the reconstruction. 
385   Int_t nRecon = 0;
386   // Loop over the detectors 
387   for (Int_t i = 1; i <= 3; i++) {
388     AliFMDSubDetector* sub = 0;
389     switch (i) {
390     case 1: sub = fFMD->GetFMD1(); break;
391     case 2: sub = fFMD->GetFMD2(); break;
392     case 3: sub = fFMD->GetFMD3(); break;
393     }
394     if (!sub) continue;
395         
396     // Loop over the rings in the detector
397     for (Int_t j = 0; j < 2; j++) {
398       Float_t     rZ = 0;
399       AliFMDRing* r  = 0;
400       switch (j) {
401       case 0: r  = sub->GetInner(); rZ = sub->GetInnerZ(); break;
402       case 1: r  = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
403       }
404       if (!r) continue;
405       
406       // Calculate low/high theta and eta 
407       // FIXME: Is this right? 
408       Float_t realZ    = fCurrentVertex + rZ;
409       Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
410       Float_t thetaIn  = TMath::ATan2(r->GetLowR(), realZ);
411       Float_t etaOut   = - TMath::Log(TMath::Tan(thetaOut / 2));
412       Float_t etaIn    = - TMath::Log(TMath::Tan(thetaIn / 2));
413       if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
414         Float_t tmp = etaIn;
415         etaIn       = etaOut;
416         etaOut      = tmp;
417       }
418
419       //-------------------------------------------------------------
420       //
421       // Here starts poisson method 
422       //
423       // Calculate eta step per strip, number of eta steps, number of
424       // phi steps, and check the sign of the eta increment 
425       Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
426       Int_t   nEta     = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta); 
427       Int_t   nPhi     = Int_t(360. / fDeltaPhi);
428       Float_t sign     = TMath::Sign(Float_t(1.), etaIn);
429
430       AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
431                         sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
432
433       // Loop over relevant phi values 
434       for (Int_t p = 0; p < nPhi; p++) {
435         Float_t  minPhi    = p * fDeltaPhi;
436         Float_t  maxPhi    = minPhi + fDeltaPhi;
437         UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
438         UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
439         
440         AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
441                           minPhi, maxPhi, minSector, maxSector));
442         // Loop over relevant eta values 
443         for (Int_t e = nEta; e >= 0; --e) {
444           Float_t  maxEta   = etaIn  - sign * e * fDeltaEta;
445           Float_t  minEta   = maxEta - sign * fDeltaEta;
446           if (sign > 0)  minEta = TMath::Max(minEta, etaOut);
447           else           minEta = TMath::Min(minEta, etaOut);
448           Float_t  theta1   = 2 * TMath::ATan(TMath::Exp(-minEta));
449           Float_t  theta2   = 2 * TMath::ATan(TMath::Exp(-maxEta));
450           Float_t  minR     = TMath::Abs(realZ * TMath::Tan(theta2));
451           Float_t  maxR     = TMath::Abs(realZ * TMath::Tan(theta1));
452           UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
453           UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
454
455           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
456                             "    [radii %f, %f, thetas %f, %f, sign %d]", 
457                             minEta, maxEta, minStrip, maxStrip, 
458                             minR, maxR, theta1, theta2, sign));
459
460           // Count number of empty strips
461           Int_t   emptyStrips = 0;
462           for (Int_t sector = minSector; sector < maxSector; sector++) 
463             for (Int_t strip = minStrip; strip < maxStrip; strip++) emptyStrips++;
464           // if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip) 
465           //     < fThreshold) emptyStrips++;
466           
467           // The total number of strips 
468           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
469           
470           // Log ratio of empty to total number of strips 
471           AliDebug(10, Form("Lambda= %d / %d = %f", 
472                             emptyStrips, nTotal, 
473                             Float_t(emptyStrips) / nTotal));
474           
475           Double_t lambda = (emptyStrips > 0 ? 
476                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
477                              1);
478
479           // The reconstructed number of particles is then given by 
480           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
481             
482           // Add a AliFMDParticles to the reconstruction tree. 
483           new((*fParticles)[nRecon])   
484             AliFMDParticles(sub->GetId(), r->GetId(),
485                             minSector, maxSector, minStrip, maxStrip,
486                             minEta, maxEta, minPhi, maxPhi,
487                             reconstructed, AliFMDParticles::kPoission);
488           nRecon++;
489         } // phi 
490       } // eta
491     } // ring 
492   } // detector 
493 }
494
495  
496 //____________________________________________________________________
497 void 
498 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/, 
499                              AliESD* /*esd*/) const
500 {
501   // nothing to be done
502
503 }
504
505 //____________________________________________________________________
506 //
507 // EOF
508 //