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