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