]> git.uio.no Git - u/mrichter/AliRoot.git/blob - FMD/AliFMDReconstructor.cxx
Fixes for reading zero-suppressed data. These should be propagated to
[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 /* $Id$ */
16 /** @file    AliFMDReconstructor.cxx
17     @author  Christian Holm Christensen <cholm@nbi.dk>
18     @date    Mon Mar 27 12:47:09 2006
19     @brief   FMD reconstruction 
20 */
21 //____________________________________________________________________
22 //
23 // This is a class that constructs AliFMDRecPoint objects from of Digits
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).   The rec-points are made via the naiive
27 // method. 
28 //
29 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
30 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
31 //
32 //
33 //____________________________________________________________________
34
35 // #include <AliLog.h>                        // ALILOG_H
36 // #include <AliRun.h>                        // ALIRUN_H
37 #include "AliFMDDebug.h"
38 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
39 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
40 #include "AliFMDAltroMapping.h"            // ALIFMDALTROMAPPING_H
41 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
42 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
43 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
44 #include "AliFMDRecPoint.h"                // ALIFMDMULTNAIIVE_H
45 #include "AliESDEvent.h"                   // ALIESDEVENT_H
46 #include "AliESDVertex.h"                  // ALIESDVERTEX_H
47 #include <AliESDFMD.h>                     // ALIESDFMD_H
48 #include <TMath.h>
49 #include <TH1.h>
50 #include <TH2.h>
51 #include <TFile.h>
52 #include <climits>
53
54 class AliRawReader;
55
56 //____________________________________________________________________
57 ClassImp(AliFMDReconstructor)
58 #if 0
59   ; // This is here to keep Emacs for indenting the next line
60 #endif
61
62 //____________________________________________________________________
63 AliFMDReconstructor::AliFMDReconstructor() 
64   : AliReconstructor(),
65     fMult(0x0), 
66     fNMult(0),
67     fTreeR(0x0),
68     fCurrentVertex(0),
69     fESDObj(0x0),
70     fNoiseFactor(0),
71     fAngleCorrect(kTRUE),
72     fVertexType(kNoVertex),
73     fESD(0x0),
74     fDiagnostics(kTRUE),
75     fDiagStep1(0), 
76     fDiagStep2(0),
77     fDiagStep3(0),
78     fDiagStep4(0),
79     fDiagAll(0)
80 {
81   // Make a new FMD reconstructor object - default CTOR.  
82   SetNoiseFactor();
83   SetAngleCorrect();
84   if (AliDebugLevel() > 0) fDiagnostics = kTRUE;
85 }
86   
87
88 //____________________________________________________________________
89 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
90   : AliReconstructor(), 
91     fMult(other.fMult),
92     fNMult(other.fNMult),
93     fTreeR(other.fTreeR),
94     fCurrentVertex(other.fCurrentVertex),
95     fESDObj(other.fESDObj),
96     fNoiseFactor(other.fNoiseFactor),
97     fAngleCorrect(other.fAngleCorrect),
98     fVertexType(other.fVertexType),
99     fESD(other.fESD),
100     fDiagnostics(other.fDiagnostics),
101     fDiagStep1(other.fDiagStep1), 
102     fDiagStep2(other.fDiagStep2),
103     fDiagStep3(other.fDiagStep3),
104     fDiagStep4(other.fDiagStep4),
105     fDiagAll(other.fDiagAll) 
106 {
107   // Copy constructor 
108 }
109   
110
111 //____________________________________________________________________
112 AliFMDReconstructor&
113 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
114 {
115   // Assignment operator
116   fMult          = other.fMult;
117   fNMult         = other.fNMult;
118   fTreeR         = other.fTreeR;
119   fCurrentVertex = other.fCurrentVertex;
120   fESDObj        = other.fESDObj;
121   fNoiseFactor   = other.fNoiseFactor;
122   fAngleCorrect  = other.fAngleCorrect;
123   fVertexType    = other.fVertexType;
124   fESD           = other.fESD;
125   fDiagnostics   = other.fDiagnostics;
126   fDiagStep1     = other.fDiagStep1;
127   fDiagStep2     = other.fDiagStep2;
128   fDiagStep3     = other.fDiagStep3;
129   fDiagStep4     = other.fDiagStep4;
130   fDiagAll       = other.fDiagAll;
131   return *this;
132 }
133
134 //____________________________________________________________________
135 AliFMDReconstructor::~AliFMDReconstructor() 
136 {
137   // Destructor 
138   if (fMult)   fMult->Delete();
139   if (fMult)   delete fMult;
140   if (fESDObj) delete fESDObj;
141 }
142
143 //____________________________________________________________________
144 void 
145 AliFMDReconstructor::Init() 
146 {
147   // Initialize the reconstructor 
148
149   // Initialize the geometry 
150   AliFMDGeometry* geom = AliFMDGeometry::Instance();
151   geom->Init();
152   geom->InitTransformations();
153
154   // Initialize the parameters
155   AliFMDParameters* param = AliFMDParameters::Instance();
156   param->Init();
157   
158   // Current vertex position
159   fCurrentVertex = 0;
160   // Create array of reconstructed strip multiplicities 
161   fMult = new TClonesArray("AliFMDRecPoint", 51200);
162   // Create ESD output object 
163   fESDObj = new AliESDFMD;
164   
165   // Check if we need diagnostics histograms 
166   if (!fDiagnostics) return;
167   AliInfo("Making diagnostics histograms");
168   fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
169                         1024, -.5, 1023.5, 1024, -.5, 1023.5);
170   fDiagStep1->SetDirectory(0);
171   fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
172   fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)", 
173                                         fNoiseFactor));
174   fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
175                         1024, -.5, 1023.5, 100, 0, 2);
176   fDiagStep2->SetDirectory(0);
177   fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
178   fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
179   fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
180                         100, 0., 2., 100, 0., 2.);
181   fDiagStep3->SetDirectory(0);
182   fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
183   fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
184   fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced", 
185                         100, 0., 2., 100, -.1, 19.9);
186   fDiagStep4->SetDirectory(0);
187   fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
188   fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
189   fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced", 
190                          1024, -.5, 1023.5, 100, -.1, 19.9);
191   fDiagAll->SetDirectory(0);
192   fDiagAll->GetXaxis()->SetTitle("ADC (read)");
193   fDiagAll->GetYaxis()->SetTitle("Multiplicity");
194 }
195
196 //____________________________________________________________________
197 void 
198 AliFMDReconstructor::ConvertDigits(AliRawReader* reader, 
199                                    TTree* digitsTree) const
200 {
201   // Convert Raw digits to AliFMDDigit's in a tree 
202   AliFMDDebug(2, ("Reading raw data into digits tree"));
203   AliFMDRawReader rawRead(reader, digitsTree);
204   // rawRead.SetSampleRate(fFMD->GetSampleRate());
205   rawRead.Exec();
206   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
207   for (size_t i = 1; i <= 3; i++) { 
208     fZS[i]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
209     fZSFactor[i] = rawRead.NoiseFactor(map->Detector2DDL(i));
210   }
211 }
212
213 //____________________________________________________________________
214 void 
215 AliFMDReconstructor::GetVertex() const
216 {
217   // Return the vertex to use. 
218   // This is obtained from the ESD object. 
219   // If not found, a warning is issued.
220   fVertexType    = kNoVertex;
221   fCurrentVertex = 0;
222   if (fESD) {
223     const AliESDVertex* vertex = fESD->GetVertex();
224     if (vertex) {
225       AliFMDDebug(2, ("Got vertex from ESD: %f", vertex->GetZv()));
226       fCurrentVertex = vertex->GetZv();
227       fVertexType    = kESDVertex;
228       return;
229     }
230   }
231   AliWarning("Didn't get any vertex from ESD or generator");
232 }
233   
234
235 //____________________________________________________________________
236 void 
237 AliFMDReconstructor::Reconstruct(AliRawReader* /*reader*/, TTree*) const
238 {
239   // Reconstruct directly from raw data (no intermediate output on
240   // digit tree or rec point tree).  
241   // Parameters: 
242   //   reader   Raw event reader 
243   //   ctree    Not used. 
244   AliError("Method is not used");
245 #if 0
246   TClonesArray*   array = new TClonesArray("AliFMDDigit");
247   AliFMDRawReader rawRead(reader, 0);
248   rawRead.ReadAdcs(array);
249   ProcessDigits(array);
250   array->Delete();
251   delete array;
252 #endif
253 }
254
255 //____________________________________________________________________
256 void 
257 AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
258                                  TTree* clusterTree) const 
259 {
260   // Reconstruct event from digits in tree 
261   // Get the FMD branch holding the digits. 
262   // FIXME: The vertex may not be known yet, so we may have to move
263   // some of this to FillESD. 
264   AliFMDDebug(2, ("Reconstructing from digits in a tree"));
265   GetVertex();
266
267   TBranch *digitBranch = digitsTree->GetBranch("FMD");
268   if (!digitBranch) {
269     Error("Exec", "No digit branch for the FMD found");
270     return;
271   }
272   TClonesArray* digits = new TClonesArray("AliFMDDigit");
273   digitBranch->SetAddress(&digits);
274
275   if (fMult)   fMult->Clear();
276   if (fESDObj) fESDObj->Clear();
277   
278   fNMult = 0;
279   fTreeR = clusterTree;
280   fTreeR->Branch("FMD", &fMult);
281   
282   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
283   digitBranch->GetEntry(0);
284   
285   AliFMDDebug(5, ("Processing digits"));
286   ProcessDigits(digits);
287
288   Int_t written = clusterTree->Fill();
289   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
290   digits->Delete();
291   delete digits;
292 }
293  
294
295 //____________________________________________________________________
296 void
297 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
298 {
299   // For each digit, find the pseudo rapdity, azimuthal angle, and
300   // number of corrected ADC counts, and pass it on to the algorithms
301   // used. 
302   Int_t nDigits = digits->GetEntries();
303   AliFMDDebug(1, ("Got %d digits", nDigits));
304   fESDObj->SetNoiseFactor(fNoiseFactor);
305   fESDObj->SetAngleCorrected(fAngleCorrect);
306   for (Int_t i = 0; i < nDigits; i++) {
307     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
308     AliFMDParameters* param  = AliFMDParameters::Instance();
309     // Check that the strip is not marked as dead 
310     if (param->IsDead(digit->Detector(), digit->Ring(), 
311                       digit->Sector(), digit->Strip())) {
312       AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", digit->Detector(), 
313                         digit->Ring(), digit->Sector(), digit->Strip()));
314       continue;
315     }
316
317     // digit->Print();
318     // Get eta and phi 
319     Float_t eta, phi;
320     PhysicalCoordinates(digit, eta, phi);
321     
322     // Substract pedestal. 
323     UShort_t counts   = SubtractPedestal(digit);
324     if(counts == USHRT_MAX) continue;
325     
326     // Gain match digits. 
327     Double_t edep     = Adc2Energy(digit, eta, counts);
328     // Get rid of nonsense energy
329     if(edep < 0)  continue;
330         
331     // Make rough multiplicity 
332     Double_t mult     = Energy2Multiplicity(digit, edep);
333     // Get rid of nonsense mult
334     if(mult < 0)  continue; 
335     AliFMDDebug(5, ("FMD%d%c[%2d,%3d]: "
336                       "ADC: %d, Counts: %d, Energy: %f, Mult: %f", 
337                       digit->Detector(), digit->Ring(), digit->Sector(),
338                       digit->Strip(), digit->Counts(), counts, edep, mult));
339     
340     // Create a `RecPoint' on the output branch. 
341     if (fMult) {
342       AliFMDRecPoint* m = 
343         new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(), 
344                                               digit->Ring(), 
345                                               digit->Sector(),
346                                               digit->Strip(),
347                                               eta, phi, 
348                                               edep, mult);
349       (void)m; // Suppress warnings about unused variables. 
350       fNMult++;
351     }
352     
353     fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(), 
354                              digit->Sector(),  digit->Strip(), mult);
355     fESDObj->SetEta(digit->Detector(), digit->Ring(), 
356                     digit->Sector(),  digit->Strip(), eta);
357
358     if (fDiagAll) fDiagAll->Fill(digit->Counts(), mult);  
359   }
360 }
361
362 //____________________________________________________________________
363 UShort_t
364 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
365 {
366   // Member function to subtract the pedestal from a digit
367   // This implementation does nothing, but a derived class could over
368   // load this to subtract a pedestal that was given in a database or
369   // something like that. 
370
371   AliFMDParameters* param  = AliFMDParameters::Instance();
372   Bool_t            zs     = fZS[digit->Detector()-1];
373   UShort_t          fac    = fZSFactor[digit->Detector()-1];
374   Float_t           ped    = (zs ? 0 : 
375                               param->GetPedestal(digit->Detector(), 
376                                                  digit->Ring(), 
377                                                  digit->Sector(), 
378                                                  digit->Strip()));
379   Float_t           noise  = param->GetPedestalWidth(digit->Detector(), 
380                                                      digit->Ring(), 
381                                                      digit->Sector(), 
382                                                      digit->Strip());
383   if(ped < 0 || noise < 0) { 
384     AliWarning(Form("Invalid pedestal (%f) or noise (%f) for %s", 
385                     ped, noise, digit->GetName()));
386     return USHRT_MAX;
387   }
388
389   AliFMDDebug(5, ("Subtracting pedestal %f from signal %d", 
390                   ped, digit->Counts()));
391   // if (digit->Count3() > 0)      adc = digit->Count3();
392   // else if (digit->Count2() > 0) adc = digit->Count2();
393   // else                          adc = digit->Count1();
394   Int_t adc    = digit->Counts();
395   Int_t counts = adc + Int_t(zs ? fac * noise : - ped);
396   counts       = TMath::Max(Int_t(counts), 0);
397   if (counts < noise * fNoiseFactor) counts = 0;
398   if (counts > 0) AliFMDDebug(15, ("Got a hit strip"));
399   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
400   
401   return  UShort_t(counts);
402 }
403
404 //____________________________________________________________________
405 Float_t
406 AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit, 
407                                 Float_t      eta, 
408                                 UShort_t     count) const
409 {
410   // Converts number of ADC counts to energy deposited. 
411   // Note, that this member function can be overloaded by derived
412   // classes to do strip-specific look-ups in databases or the like,
413   // to find the proper gain for a strip. 
414   // 
415   // In the first simple version, we calculate the energy deposited as 
416   // 
417   //    EnergyDeposited = cos(theta) * gain * count
418   // 
419   // where 
420   // 
421   //           Pre_amp_MIP_Range
422   //    gain = ----------------- * Energy_deposited_per_MIP
423   //           ADC_channel_size    
424   // 
425   // is constant and the same for all strips.
426   //
427   // For the production we use the conversion measured in the NBI lab.
428   // The total conversion is then:
429   // gain = ADC / DAC
430   // => energy = EdepMip*count / gain*DACPerADC
431   //
432   //
433   if (count <= 0) return 0;
434   AliFMDParameters* param = AliFMDParameters::Instance();
435   Float_t           gain  = param->GetPulseGain(digit->Detector(), 
436                                                 digit->Ring(), 
437                                                 digit->Sector(), 
438                                                 digit->Strip());
439   // 'Tagging' bad gains as bad energy
440   if (gain < 0) { 
441     AliWarning(Form("Invalid gain (%f) for %s", gain, digit->GetName()));
442     return -1;
443   }
444   AliFMDDebug(5, ("Converting counts %d to energy via factor %f and DAC2MIP %f", 
445                   count, gain,param->GetDACPerMIP()));
446
447   Double_t edep  = (count * param->GetEdepMip()) / (gain * param->GetDACPerMIP());
448   if (fDiagStep2) fDiagStep2->Fill(count, edep);  
449   if (fAngleCorrect) {
450     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
451     Double_t corr  = TMath::Abs(TMath::Cos(theta));
452     Double_t cedep = corr * edep;
453     AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
454                       edep, corr, cedep, eta, theta));
455     if (fDiagStep3) fDiagStep3->Fill(edep, cedep);  
456     edep = cedep;
457   }
458   return edep;
459 }
460
461 //____________________________________________________________________
462 Float_t
463 AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */, 
464                                          Float_t      edep) const
465 {
466   // Converts an energy signal to number of particles. 
467   // Note, that this member function can be overloaded by derived
468   // classes to do strip-specific look-ups in databases or the like,
469   // to find the proper gain for a strip. 
470   // 
471   // In this simple version, we calculate the multiplicity as 
472   // 
473   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
474   // 
475   // where 
476   //
477   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness 
478   // 
479   // is constant and the same for all strips 
480   AliFMDParameters* param   = AliFMDParameters::Instance();
481   Double_t          edepMIP = param->GetEdepMip();
482   Float_t           mult    = edep / edepMIP;
483 #if 0
484   if (edep > 0) 
485     AliFMDDebug(15, ("Translating energy %f to multiplicity via "
486                      "divider %f->%f", edep, edepMIP, mult));
487 #endif
488   if (fDiagStep4) fDiagStep4->Fill(edep, mult);  
489   return mult;
490 }
491
492 //____________________________________________________________________
493 void
494 AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit, 
495                                          Float_t& eta, 
496                                          Float_t& phi) const
497 {
498   // Get the eta and phi of a digit 
499   // 
500   // Get geometry. 
501   AliFMDGeometry* geom = AliFMDGeometry::Instance();
502   Double_t x, y, z, r, theta;
503   geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(), 
504                     digit->Strip(), x, y, z);
505   // Correct for vertex offset. 
506   z     += fCurrentVertex;
507   phi   =  TMath::ATan2(y, x);
508   r     =  TMath::Sqrt(y * y + x * x);
509   theta =  TMath::ATan2(r, z);
510   eta   = -TMath::Log(TMath::Tan(theta / 2));
511 }
512
513       
514
515 //____________________________________________________________________
516 void 
517 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
518                              TTree*  /* clusterTree */,
519                              AliESDEvent* esd) const
520 {
521   // nothing to be done
522   // FIXME: The vertex may not be known when Reconstruct is executed,
523   // so we may have to move some of that member function here. 
524   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
525   // fESDObj->Print();
526
527   if (esd) { 
528     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
529     esd->SetFMDData(fESDObj);
530   }
531
532   if (!fDiagnostics || !esd) return;
533   static bool first = true;
534   // This is most likely NOT the event number you'd like to use. It
535   // has nothing to do with the 'real' event number. 
536   // - That's OK.  We just use it for the name of the directory -
537   // nothing else.  Christian
538   Int_t evno = esd->GetEventNumberInFile(); 
539   AliFMDDebug(1, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
540   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
541   first = false;
542   f.cd(); 
543   TDirectory* d = f.mkdir(Form("%03d", evno), 
544                           Form("Diagnostics histograms for event # %d", evno));
545   d->cd();
546   if (fDiagStep1) fDiagStep1->Write();
547   if (fDiagStep2) fDiagStep2->Write();
548   if (fDiagStep3) fDiagStep3->Write();
549   if (fDiagStep4) fDiagStep4->Write();
550   if (fDiagAll)   fDiagAll->Write();
551   d->Write();
552   f.Write();
553   f.Close();
554
555   if (fDiagStep1) fDiagStep1->Reset();
556   if (fDiagStep2) fDiagStep2->Reset();
557   if (fDiagStep3) fDiagStep3->Reset();
558   if (fDiagStep4) fDiagStep4->Reset();
559   if (fDiagAll)   fDiagAll->Reset();
560 }
561
562 //____________________________________________________________________
563 void
564 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree, 
565                              AliESDEvent* esd) const
566 {
567   TTree* dummy = 0;
568   FillESD(dummy, clusterTree, esd);
569 }
570
571 //____________________________________________________________________
572 //
573 // EOF
574 //