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