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