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