0468f0b0e7dbc9253a517cbb8266f7e4095c9037
[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"
47 #include "AliFMDDigit.h"
48 #include "AliFMDParticles.h"
49 #include "AliFMDReconstructor.h"
50 #include "AliAltroBuffer.h"
51 #include "AliLog.h"
52 #include "AliRun.h"
53 #include "AliRunLoader.h"
54 #include "AliLoader.h"
55 #include "AliHeader.h"
56 #include "AliGenEventHeader.h"
57 #include "AliFMDRawStream.h"
58 #include "AliRawReader.h"
59
60 //____________________________________________________________________
61 ClassImp(AliFMDReconstructor);
62
63 //____________________________________________________________________
64 AliFMDReconstructor::AliFMDReconstructor() 
65   : AliReconstructor(),
66     fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
67     fDeltaEta(0), 
68     fDeltaPhi(0), 
69     fThreshold(0),
70     fPedestal(0), 
71     fPedestalWidth(0)
72 {
73   // Make a new FMD reconstructor object - default CTOR.
74   SetDeltaEta();
75   SetDeltaPhi();
76   SetThreshold();
77   SetPedestal();
78
79   fParticles = new TClonesArray("AliFMDParticles", 1000);
80   fFMDLoader = 0;
81   fRunLoader = 0;
82   fFMD       = 0;
83 }
84   
85
86 //____________________________________________________________________
87 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
88   : AliReconstructor(),
89     fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
90     fDeltaEta(0), 
91     fDeltaPhi(0), 
92     fThreshold(0),
93     fPedestal(0), 
94     fPedestalWidth(0)
95 {
96   // Make a new FMD reconstructor object - default CTOR.
97   SetDeltaEta(other.fDeltaEta);
98   SetDeltaPhi(other.fDeltaPhi);
99   SetThreshold(other.fThreshold);
100   SetPedestal(other.fPedestal, other.fPedestalWidth);
101
102   // fParticles = new TClonesArray("AliFMDParticles", 1000);
103   fFMDLoader = other.fFMDLoader;
104   fRunLoader = other.fRunLoader;
105   fFMD       = other.fFMD;
106 }
107   
108
109 //____________________________________________________________________
110 AliFMDReconstructor&
111 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
112 {
113   // Make a new FMD reconstructor object - default CTOR.
114   SetDeltaEta(other.fDeltaEta);
115   SetDeltaPhi(other.fDeltaPhi);
116   SetThreshold(other.fThreshold);
117   SetPedestal(other.fPedestal, other.fPedestalWidth);
118
119   // fParticles = new TClonesArray("AliFMDParticles", 1000);
120   fFMDLoader = other.fFMDLoader;
121   fRunLoader = other.fRunLoader;
122   fFMD       = other.fFMD;
123
124   return *this;
125 }
126   
127 //____________________________________________________________________
128 void 
129 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width) 
130 {
131   // Set the pedestal, and pedestal width 
132   fPedestal      = mean;
133   fPedestalWidth = width;
134 }
135
136 //____________________________________________________________________
137 void 
138 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader, 
139                                  AliRawReader* rawReader) const
140
141   // Collects all digits in the same active volume into number of
142   // particles
143   //
144   // Reconstruct number of particles in given group of pads for given
145   // FMDvolume determined by numberOfVolume,
146   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
147   // numberOgMaxRing 
148   //
149   // The reconstruction method is choosen based on the number of empty
150   // strips. 
151   fParticles->Clear();
152   if (!runLoader) {
153     Error("Exec","Run Loader loader is NULL - Session not opened");
154     return;
155   }
156   fRunLoader = runLoader;
157   fFMDLoader = runLoader->GetLoader("FMDLoader");
158   if (!fFMDLoader) 
159     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
160           "in specified event");
161
162   // Get the AliRun object
163   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
164   
165   // Get the AliFMD object
166   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
167   if (!fFMD) {
168     AliError("Can not get FMD from gAlice");
169     return;
170   }
171   fFMDLoader->LoadRecPoints("RECREATE");
172
173   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
174
175   TClonesArray* digits = fFMD->Digits();
176   if (rawReader) {
177     Int_t event = 0;
178     while (rawReader->NextEvent()) {
179       ProcessEvent(event, rawReader, digits);
180       event++;
181     }
182   }
183   else {
184     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries()); 
185     for(Int_t event = 0; event < nEvents; event++) 
186       ProcessEvent(event, 0, digits);
187   }
188
189
190   fFMDLoader->UnloadRecPoints();
191   fFMDLoader = 0;
192   fRunLoader = 0;
193   fFMD       = 0;
194 }
195
196 //____________________________________________________________________
197 void 
198 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
199
200   // Collects all digits in the same active volume into number of
201   // particles
202   //
203   // Reconstruct number of particles in given group of pads for given
204   // FMDvolume determined by numberOfVolume,
205   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
206   // numberOgMaxRing 
207   //
208   // The reconstruction method is choosen based on the number of empty
209   // strips. 
210   Reconstruct(runLoader, 0);
211 }
212
213
214 //____________________________________________________________________
215 void 
216 AliFMDReconstructor::ProcessEvent(Int_t event, 
217                                   AliRawReader* reader, 
218                                   TClonesArray* digits) const
219 {
220   // Process one event read from either a clones array or from a a raw
221   // data reader. 
222   fRunLoader->GetEvent(event) ;
223   //event z-vertex for correction eta-rad dependence      
224   AliHeader *header            = fRunLoader->GetHeader();
225   if (!header) 
226     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   Float_t zVertex = o.At(2);
233
234   // If the recontruction tree isn't loaded, load it
235   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
236   
237   //Make branches to hold the reconstructed particles 
238   const Int_t kBufferSize = 16000;
239   fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
240
241   // Load or recreate the digits 
242   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
243     if (!reader) {
244       Error("Exec","Error occured while loading digits. Exiting.");
245       return;
246     }
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   if (!digitBranch) {
264     if (!reader) {
265       Error("Exec", "No digit branch for the FMD found");
266       return;
267     }
268     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
269   }
270   if (!reader) digitBranch->SetAddress(&digits);
271
272   fEmptyStrips = 0;
273   fTotalStrips = 0;
274   Bool_t ok = kFALSE;
275   if      (reader) ok = ReadAdcs(reader);
276   else if (digits) ok = ReadAdcs(digits);
277   if (!ok) return;
278   
279   ReconstructFromCache(zVertex);
280
281   if (reader) {
282     digitTree->Fill();
283     fFMDLoader->WriteDigits("OVERWRITE");
284   }
285   fFMDLoader->UnloadDigits();
286   fFMDLoader->TreeR()->Reset();
287   fFMDLoader->TreeR()->Fill(); 
288   fFMDLoader->WriteRecPoints("OVERWRITE");
289 }
290
291 //____________________________________________________________________
292 Bool_t
293 AliFMDReconstructor::ReadAdcs(TClonesArray* digits) const
294 {
295   // Read the ADC values from a clones array. 
296   AliDebug(10, "Reading ADCs from Digits array");
297   // read Digits, and reconstruct the particles
298   if (!fFMDLoader->TreeD()->GetEvent(0)) return kFALSE;
299
300   // Reads the digits from the array, and fills up the cache (fAdcs) 
301   fAdcs.Clear();
302   Int_t nDigits = digits->GetEntries();
303   for (Int_t digit = 0; digit < nDigits; digit++) {
304     AliFMDDigit* fmdDigit = 
305       static_cast<AliFMDDigit*>(digits->UncheckedAt(digit));    
306     
307     ProcessDigit(fmdDigit);
308   } //digit loop
309   return kTRUE;
310 }
311
312 //____________________________________________________________________
313 Bool_t
314 AliFMDReconstructor::ReadAdcs(AliRawReader* reader) const
315 {
316   // Read the ADC values from a raw data reader. 
317   AliDebug(10, "Reading ADCs from RawReader");
318   // Reads the digits from a RAW data 
319   fAdcs.Clear();
320   // reader->Reset();
321
322   if (!reader->ReadHeader()) {
323     Error("ReadAdcs", "Couldn't read header");
324     return kFALSE;
325   }
326
327   // Use AliAltroRawStream to read the ALTRO format.  No need to
328   // reinvent the wheel :-) 
329   AliFMDRawStream input(reader);
330   // Select FMD DDL's 
331   reader->Select(AliFMD::kBaseDDL >> 8);
332   
333   Int_t    oldDDL      = -1;
334   Int_t    count       = 0;
335   UShort_t detector    = 1; // Must be one here
336   UShort_t oldDetector = 0;
337   // Loop over data in file 
338   Bool_t   next       = kTRUE;
339
340   // local Cache 
341   TArrayI counts(10);
342   counts.Reset(-1);
343   Int_t offset = 0;
344   
345   while (next) {
346     next = input.Next();
347
348
349     count++; 
350     Int_t ddl = reader->GetDDLID();
351     if (ddl != oldDDL 
352         || input.IsNewStrip()
353         || !next) {
354       // Make a new digit, if we have some data! 
355       if (counts[0] >= 0) {
356         // Got a new strip. 
357         AliDebug(10, Form("Add a new strip: FMD%d%c[%2d,%3d] "
358                           "(current: FMD%d%c[%2d,%3d])", 
359                           oldDetector, input.PrevRing(), 
360                           input.PrevSector() , input.PrevStrip(),
361                           detector , input.Ring(), input.Sector(), 
362                           input.Strip()));
363         fFMD->AddDigit(oldDetector, 
364                        input.PrevRing(), 
365                        input.PrevSector(), 
366                        input.PrevStrip(), 
367                        counts[0], counts[1], counts[2]);
368         AliFMDDigit* digit = 
369           static_cast<AliFMDDigit*>(fFMD->Digits()->
370                                     UncheckedAt(fFMD->GetNdigits()-1));
371         ProcessDigit(digit);
372       }
373         
374       if (!next) { 
375         AliDebug(10, Form("Read %d channels for FMD%d", 
376                           count + 1, detector));
377         break;
378       }
379     
380     
381       // If we got a new DDL, it means we have a new detector. 
382       if (ddl != oldDDL) {
383         if (detector != 0) 
384           AliDebug(10, Form("Read %d channels for FMD%d", 
385                             count + 1, detector));
386         // Reset counts, and update the DDL cache 
387         count       = 0;
388         oldDDL      = ddl;
389         // Check that we're processing a FMD detector 
390         Int_t detId = reader->GetDetectorID();
391         if (detId != (AliFMD::kBaseDDL >> 8)) {
392           Error("ReadAdcs", "Detector ID %d != %d",
393                 detId, (AliFMD::kBaseDDL >> 8));
394           break;
395         }
396         // Figure out what detector we're deling with 
397         oldDetector = detector;
398         switch (ddl) {
399         case 0: detector = 1; break;
400         case 1: detector = 2; break;
401         case 2: detector = 3; break;
402         default:
403           Error("ReadAdcs", "Unknown DDL 0x%x for FMD", ddl);
404           return kFALSE;
405         }
406         AliDebug(10, Form("Reading ADCs for 0x%x  - That is FMD%d",
407                           reader->GetEquipmentId(), detector));
408       }
409       counts.Reset(-1);
410       offset = 0;
411     }
412     
413     counts[offset] = input.Count();
414     offset++;
415     
416     AliDebug(10, Form("ADC of FMD%d%c[%2d,%3d] += %d",
417                       detector, input.Ring(), input.Sector(), 
418                       input.Strip(), input.Count()));
419     oldDetector = detector;
420   }
421   return kTRUE;
422 }
423
424 //____________________________________________________________________
425 void
426 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
427 {
428   // Process a digit.  Derived classes can overload this member
429   // function to do stuff to the digit.  However, it should write the
430   // ADC count to the internal cache 
431   //
432   //    fAdcs(detector - 1, ring, sector, strip) = counts;
433   //
434   // In this implementation, we count the number of strips below
435   // threshold.   This we do to later choose what kind of
436   // reconstruction algorithm we'd like to use. 
437   // 
438   UShort_t detector = digit->Detector();
439   Char_t   ring     = digit->Ring();
440   UShort_t sector   = digit->Sector();
441   UShort_t strip    = digit->Strip();    
442   
443   UShort_t counts  = SubtractPedestal(digit);
444   
445   fAdcs(detector - 1, ring, sector, strip) = counts;
446   if (counts < fThreshold) fEmptyStrips++;
447   fTotalStrips++;
448 }
449
450 //____________________________________________________________________
451 UShort_t
452 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
453 {
454   // Member function to subtract the pedestal from a digit
455   // This implementation does nothing, but a derived class could over
456   // load this to subtract a pedestal that was given in a database or
457   // something like that. 
458
459   Int_t counts = 
460     TMath::Max(Int_t(digit->Count1() - fPedestal - 3 * fPedestalWidth), 0);
461   if (digit->Count2() >= 0) 
462     counts += 
463       TMath::Max(Int_t(digit->Count2() - fPedestal - 3 * fPedestalWidth), 0);
464   if (digit->Count3() >= 0) 
465     counts += 
466       TMath::Max(Int_t(digit->Count3() - fPedestal - 3 * fPedestalWidth), 0);
467       
468   if (counts < 0) counts = 0;
469   return  UShort_t(counts);
470 }
471
472 //____________________________________________________________________
473 void
474 AliFMDReconstructor::ReconstructFromCache(Float_t zVertex) const
475 {
476   // Based on the information in the cache, do the reconstruction. 
477   Int_t nRecon = 0;
478   // Loop over the detectors 
479   for (Int_t i = 1; i <= 3; i++) {
480     AliFMDSubDetector* sub = 0;
481     switch (i) {
482     case 1: sub = fFMD->GetFMD1(); break;
483     case 2: sub = fFMD->GetFMD2(); break;
484     case 3: sub = fFMD->GetFMD3(); break;
485     }
486     if (!sub) continue;
487         
488     // Loop over the rings in the detector
489     for (Int_t j = 0; j < 2; j++) {
490       Float_t     rZ = 0;
491       AliFMDRing* r  = 0;
492       switch (j) {
493       case 0: r  = sub->GetInner(); rZ = sub->GetInnerZ(); break;
494       case 1: r  = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
495       }
496       if (!r) continue;
497       
498       // Calculate low/high theta and eta 
499       // FIXME: Is this right? 
500       Float_t realZ    = zVertex + rZ;
501       Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
502       Float_t thetaIn  = TMath::ATan2(r->GetLowR(), realZ);
503       Float_t etaOut   = - TMath::Log(TMath::Tan(thetaOut / 2));
504       Float_t etaIn    = - TMath::Log(TMath::Tan(thetaIn / 2));
505       if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
506         Float_t tmp = etaIn;
507         etaIn       = etaOut;
508         etaOut      = tmp;
509       }
510
511       //-------------------------------------------------------------
512       //
513       // Here starts poisson method 
514       //
515       // Calculate eta step per strip, number of eta steps, number of
516       // phi steps, and check the sign of the eta increment 
517       Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
518       Int_t   nEta     = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta); 
519       Int_t   nPhi     = Int_t(360. / fDeltaPhi);
520       Float_t sign     = TMath::Sign(Float_t(1.), etaIn);
521
522       AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
523                         sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
524
525       // Loop over relevant phi values 
526       for (Int_t p = 0; p < nPhi; p++) {
527         Float_t  minPhi    = p * fDeltaPhi;
528         Float_t  maxPhi    = minPhi + fDeltaPhi;
529         UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
530         UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
531         
532         AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
533                           minPhi, maxPhi, minSector, maxSector));
534         // Loop over relevant eta values 
535         for (Int_t e = nEta; e >= 0; --e) {
536           Float_t  maxEta   = etaIn  - sign * e * fDeltaEta;
537           Float_t  minEta   = maxEta - sign * fDeltaEta;
538           if (sign > 0)  minEta = TMath::Max(minEta, etaOut);
539           else           minEta = TMath::Min(minEta, etaOut);
540           Float_t  theta1   = 2 * TMath::ATan(TMath::Exp(-minEta));
541           Float_t  theta2   = 2 * TMath::ATan(TMath::Exp(-maxEta));
542           Float_t  minR     = TMath::Abs(realZ * TMath::Tan(theta2));
543           Float_t  maxR     = TMath::Abs(realZ * TMath::Tan(theta1));
544           UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
545           UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
546
547           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
548                             "    [radii %f, %f, thetas %f, %f, sign %d]", 
549                             minEta, maxEta, minStrip, maxStrip, 
550                             minR, maxR, theta1, theta2, sign));
551
552           // Count number of empty strips
553           Int_t   emptyStrips = 0;
554           for (Int_t sector = minSector; sector < maxSector; sector++) 
555             for (Int_t strip = minStrip; strip < maxStrip; strip++) 
556               if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip) 
557                   < fThreshold) emptyStrips++;
558           
559           // The total number of strips 
560           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
561           
562           // Log ratio of empty to total number of strips 
563           AliDebug(10, Form("Lambda= %d / %d = %f", 
564                             emptyStrips, nTotal, 
565                             Float_t(emptyStrips) / nTotal));
566           
567           Double_t lambda = (emptyStrips > 0 ? 
568                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
569                              1);
570
571           // The reconstructed number of particles is then given by 
572           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
573             
574           // Add a AliFMDParticles to the reconstruction tree. 
575           new((*fParticles)[nRecon])   
576             AliFMDParticles(sub->GetId(), r->GetId(),
577                             minSector, maxSector, minStrip, maxStrip,
578                             minEta, maxEta, minPhi, maxPhi,
579                             reconstructed, AliFMDParticles::kPoission);
580           nRecon++;
581         } // phi 
582       } // eta
583     } // ring 
584   } // detector 
585 }
586
587  
588 //____________________________________________________________________
589 void 
590 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/, 
591                              AliESD* /*esd*/) const
592 {
593   // nothing to be done
594
595 }
596
597 //____________________________________________________________________
598 //
599 // EOF
600 //