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