[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id$ */
/** @file    AliFMDReconstructor.cxx
    @author  Christian Holm Christensen <cholm@nbi.dk>
    @date    Mon Mar 27 12:47:09 2006
    @brief   FMD reconstruction 
*/
//____________________________________________________________________
//
// This is a class that constructs AliFMDRecPoint objects from of Digits
// This class reads either digits from a TClonesArray or raw data from 
// a DDL file (or similar), and stores the read ADC counts in an
// internal cache (fAdcs).   The rec-points are made via the naiive
// method. 
//
//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
//  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
//
//
//____________________________________________________________________

// #include <AliLog.h>                        // ALILOG_H
// #include <AliRun.h>                        // ALIRUN_H
#include "AliFMDDebug.h"
#include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
#include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
#include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
#include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
#include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
#include "AliFMDRecPoint.h"	   	   // ALIFMDMULTNAIIVE_H
#include "AliESDEvent.h"		   // ALIESDEVENT_H
#include "AliESDVertex.h"		   // ALIESDVERTEX_H
#include <AliESDFMD.h>			   // ALIESDFMD_H
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
class AliRawReader;

//____________________________________________________________________
ClassImp(AliFMDReconstructor)
#if 0
  ; // This is here to keep Emacs for indenting the next line
#endif

//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor() 
  : AliReconstructor(),
    fMult(0x0), 
    fNMult(0),
    fTreeR(0x0),
    fCurrentVertex(0),
    fESDObj(0x0),
    fNoiseFactor(0),
    fAngleCorrect(kTRUE),
    fVertexType(kNoVertex),
    fESD(0x0),
    fDiagnostics(kFALSE),
    fDiagStep1(0), 
    fDiagStep2(0),
    fDiagStep3(0),
    fDiagStep4(0),
    fDiagAll(0)
{
  // Make a new FMD reconstructor object - default CTOR.  
  SetNoiseFactor();
  SetAngleCorrect();
}
  

//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
  : AliReconstructor(), 
    fMult(other.fMult),
    fNMult(other.fNMult),
    fTreeR(other.fTreeR),
    fCurrentVertex(other.fCurrentVertex),
    fESDObj(other.fESDObj),
    fNoiseFactor(other.fNoiseFactor),
    fAngleCorrect(other.fAngleCorrect),
    fVertexType(other.fVertexType),
    fESD(other.fESD),
    fDiagnostics(other.fDiagnostics),
    fDiagStep1(other.fDiagStep1), 
    fDiagStep2(other.fDiagStep2),
    fDiagStep3(other.fDiagStep3),
    fDiagStep4(other.fDiagStep4),
    fDiagAll(other.fDiagAll) 
{
  // Copy constructor 
}
  

//____________________________________________________________________
AliFMDReconstructor&
AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
{
  // Assignment operator
  fMult          = other.fMult;
  fNMult         = other.fNMult;
  fTreeR         = other.fTreeR;
  fCurrentVertex = other.fCurrentVertex;
  fESDObj        = other.fESDObj;
  fNoiseFactor   = other.fNoiseFactor;
  fAngleCorrect  = other.fAngleCorrect;
  fVertexType    = other.fVertexType;
  fESD           = other.fESD;
  fDiagnostics   = other.fDiagnostics;
  fDiagStep1     = other.fDiagStep1;
  fDiagStep2     = other.fDiagStep2;
  fDiagStep3     = other.fDiagStep3;
  fDiagStep4     = other.fDiagStep4;
  fDiagAll       = other.fDiagAll;
  return *this;
}

//____________________________________________________________________
AliFMDReconstructor::~AliFMDReconstructor() 
{
  // Destructor 
  if (fMult)   fMult->Delete();
  if (fMult)   delete fMult;
  if (fESDObj) delete fESDObj;
}

//____________________________________________________________________
void 
AliFMDReconstructor::Init() 
{
  // Initialize the reconstructor 

  // Initialize the geometry 
  AliFMDGeometry* geom = AliFMDGeometry::Instance();
  geom->Init();
  geom->InitTransformations();

  // Initialize the parameters
  AliFMDParameters* param = AliFMDParameters::Instance();
  param->Init();
  
  // Current vertex position
  fCurrentVertex = 0;
  // Create array of reconstructed strip multiplicities 
  fMult = new TClonesArray("AliFMDRecPoint", 51200);
  // Create ESD output object 
  fESDObj = new AliESDFMD;
  
  // Check if we need diagnostics histograms 
  if (!fDiagnostics) return;
  fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
			1024, -.5, 1023.5, 1024, -.5, 1023.5);
  fDiagStep1->SetDirectory(0);
  fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
  fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)", 
					fNoiseFactor));
  fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
			1024, -.5, 1023.5, 100, 0, 2);
  fDiagStep2->SetDirectory(0);
  fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
  fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
  fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
			100, 0., 2., 100, 0., 2.);
  fDiagStep3->SetDirectory(0);
  fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
  fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
  fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced", 
			100, 0., 2., 100, -.1, 19.9);
  fDiagStep4->SetDirectory(0);
  fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
  fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
  fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced", 
			 1024, -.5, 1023.5, 100, -.1, 19.9);
  fDiagAll->SetDirectory(0);
  fDiagAll->GetXaxis()->SetTitle("ADC (read)");
  fDiagAll->GetYaxis()->SetTitle("Multiplicity");
}

//____________________________________________________________________
void 
AliFMDReconstructor::ConvertDigits(AliRawReader* reader, 
				   TTree* digitsTree) const
{
  // Convert Raw digits to AliFMDDigit's in a tree 
  AliFMDDebug(2, ("Reading raw data into digits tree"));
  AliFMDRawReader rawRead(reader, digitsTree);
  // rawRead.SetSampleRate(fFMD->GetSampleRate());
  rawRead.Exec();
}

//____________________________________________________________________
void 
AliFMDReconstructor::GetVertex() const
{
  // Return the vertex to use. 
  // This is obtained from the ESD object. 
  // If not found, a warning is issued.
  fVertexType    = kNoVertex;
  fCurrentVertex = 0;
  if (fESD) {
    const AliESDVertex* vertex = fESD->GetVertex();
    if (vertex) {
      AliFMDDebug(2, ("Got vertex from ESD: %f", vertex->GetZv()));
      fCurrentVertex = vertex->GetZv();
      fVertexType    = kESDVertex;
      return;
    }
  }
  AliWarning("Didn't get any vertex from ESD or generator");
}
  

//____________________________________________________________________
void 
AliFMDReconstructor::Reconstruct(AliRawReader* /*reader*/, TTree*) const
{
  // Reconstruct directly from raw data (no intermediate output on
  // digit tree or rec point tree).  
  // Parameters: 
  //   reader	Raw event reader 
  //   ctree    Not used. 
  AliError("Method is not used");
#if 0
  TClonesArray*   array = new TClonesArray("AliFMDDigit");
  AliFMDRawReader rawRead(reader, 0);
  rawRead.ReadAdcs(array);
  ProcessDigits(array);
  array->Delete();
  delete array;
#endif
}

//____________________________________________________________________
void 
AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
				 TTree* clusterTree) const 
{
  // Reconstruct event from digits in tree 
  // Get the FMD branch holding the digits. 
  // FIXME: The vertex may not be known yet, so we may have to move
  // some of this to FillESD. 
  AliFMDDebug(2, ("Reconstructing from digits in a tree"));
  GetVertex();

  TBranch *digitBranch = digitsTree->GetBranch("FMD");
  if (!digitBranch) {
    Error("Exec", "No digit branch for the FMD found");
    return;
  }
  TClonesArray* digits = new TClonesArray("AliFMDDigit");
  digitBranch->SetAddress(&digits);

  if (fMult)   fMult->Clear();
  if (fESDObj) fESDObj->Clear();
  
  fNMult = 0;
  fTreeR = clusterTree;
  fTreeR->Branch("FMD", &fMult);
  
  AliFMDDebug(5, ("Getting entry 0 from digit branch"));
  digitBranch->GetEntry(0);
  
  AliFMDDebug(5, ("Processing digits"));
  ProcessDigits(digits);

  Int_t written = clusterTree->Fill();
  AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
  digits->Delete();
  delete digits;
}
 

//____________________________________________________________________
void
AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
{
  // For each digit, find the pseudo rapdity, azimuthal angle, and
  // number of corrected ADC counts, and pass it on to the algorithms
  // used. 
  Int_t nDigits = digits->GetEntries();
  AliFMDDebug(1, ("Got %d digits", nDigits));
  fESDObj->SetNoiseFactor(fNoiseFactor);
  fESDObj->SetAngleCorrected(fAngleCorrect);
  for (Int_t i = 0; i < nDigits; i++) {
    AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
    AliFMDParameters* param  = AliFMDParameters::Instance();
    // Check that the strip is not marked as dead 
    if (param->IsDead(digit->Detector(), digit->Ring(), 
		      digit->Sector(), digit->Strip())) {
      AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", digit->Detector(), 
			digit->Ring(), digit->Sector(), digit->Strip()));
      continue;
    }

    // digit->Print();
    // Get eta and phi 
    Float_t eta, phi;
    PhysicalCoordinates(digit, eta, phi);
    
    // Substract pedestal. 
    UShort_t counts   = SubtractPedestal(digit);
    
    // Gain match digits. 
    Double_t edep     = Adc2Energy(digit, eta, counts);
    
    // Make rough multiplicity 
    Double_t mult     = Energy2Multiplicity(digit, edep);
    
    AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
		      "ADC: %d, Counts: %d, Energy: %f, Mult: %f", 
		      digit->Detector(), digit->Ring(), digit->Sector(),
		      digit->Strip(), digit->Counts(), counts, edep, mult));
    
    // Create a `RecPoint' on the output branch. 
    if (fMult) {
      AliFMDRecPoint* m = 
	new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(), 
					      digit->Ring(), 
					      digit->Sector(),
					      digit->Strip(),
					      eta, phi, 
					      edep, mult);
      (void)m; // Suppress warnings about unused variables. 
      fNMult++;
    }
    
    fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(), 
			     digit->Sector(),  digit->Strip(), mult);
    fESDObj->SetEta(digit->Detector(), digit->Ring(), 
		    digit->Sector(),  digit->Strip(), eta);

    if (fDiagAll) fDiagAll->Fill(digit->Counts(), mult);  
  }
}

//____________________________________________________________________
UShort_t
AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
{
  // Member function to subtract the pedestal from a digit
  // This implementation does nothing, but a derived class could over
  // load this to subtract a pedestal that was given in a database or
  // something like that. 

  Int_t             counts = 0;
  Int_t             adc    = 0;
  AliFMDParameters* param  = AliFMDParameters::Instance();
  Float_t           ped    = param->GetPedestal(digit->Detector(), 
						digit->Ring(), 
						digit->Sector(), 
						digit->Strip());
  Float_t           noise  = param->GetPedestalWidth(digit->Detector(), 
						     digit->Ring(), 
						     digit->Sector(), 
						     digit->Strip());
  AliFMDDebug(15, ("Subtracting pedestal %f from signal %d", 
		   ped, digit->Counts()));
  if (digit->Count3() > 0)      adc = digit->Count3();
  else if (digit->Count2() > 0) adc = digit->Count2();
  else                          adc = digit->Count1();
  counts = TMath::Max(Int_t(adc - ped), 0);
  if (counts < noise * fNoiseFactor) counts = 0;
  if (counts > 0) AliFMDDebug(15, ("Got a hit strip"));
  if (fDiagStep1) fDiagStep1->Fill(adc, counts);
  
  return  UShort_t(counts);
}

//____________________________________________________________________
Float_t
AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit, 
				Float_t      eta, 
				UShort_t     count) const
{
  // Converts number of ADC counts to energy deposited. 
  // Note, that this member function can be overloaded by derived
  // classes to do strip-specific look-ups in databases or the like,
  // to find the proper gain for a strip. 
  // 
  // In this simple version, we calculate the energy deposited as 
  // 
  //    EnergyDeposited = cos(theta) * gain * count
  // 
  // where 
  // 
  //           Pre_amp_MIP_Range
  //    gain = ----------------- * Energy_deposited_per_MIP
  //           ADC_channel_size    
  // 
  // is constant and the same for all strips.
  if (count <= 0) return 0;
  AliFMDParameters* param = AliFMDParameters::Instance();
  Float_t           gain  = param->GetPulseGain(digit->Detector(), 
						digit->Ring(), 
						digit->Sector(), 
						digit->Strip());
  AliFMDDebug(15, ("Converting counts %d to energy via factor %f", 
		    count, gain));

  Double_t edep  = count * gain;
  if (fDiagStep2) fDiagStep2->Fill(count, edep);  
  if (fAngleCorrect) {
    Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
    Double_t corr  = TMath::Abs(TMath::Cos(theta));
    Double_t cedep = corr * edep;
    AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
		      edep, corr, cedep, eta, theta));
    if (fDiagStep3) fDiagStep3->Fill(edep, cedep);  
    edep = cedep;
  }
  return edep;
}

//____________________________________________________________________
Float_t
AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */, 
					 Float_t      edep) const
{
  // Converts an energy signal to number of particles. 
  // Note, that this member function can be overloaded by derived
  // classes to do strip-specific look-ups in databases or the like,
  // to find the proper gain for a strip. 
  // 
  // In this simple version, we calculate the multiplicity as 
  // 
  //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
  // 
  // where 
  //
  //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness 
  // 
  // is constant and the same for all strips 
  AliFMDParameters* param   = AliFMDParameters::Instance();
  Double_t          edepMIP = param->GetEdepMip();
  Float_t           mult    = edep / edepMIP;
  if (edep > 0) 
    AliFMDDebug(15, ("Translating energy %f to multiplicity via "
		     "divider %f->%f", edep, edepMIP, mult));
  if (fDiagStep4) fDiagStep4->Fill(edep, mult);  
  return mult;
}

//____________________________________________________________________
void
AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit, 
					 Float_t& eta, 
					 Float_t& phi) const
{
  // Get the eta and phi of a digit 
  // 
  // Get geometry. 
  AliFMDGeometry* geom = AliFMDGeometry::Instance();
  Double_t x, y, z, r, theta;
  geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(), 
		    digit->Strip(), x, y, z);
  // Correct for vertex offset. 
  z     += fCurrentVertex;
  phi   =  TMath::ATan2(y, x);
  r     =  TMath::Sqrt(y * y + x * x);
  theta =  TMath::ATan2(r, z);
  eta   = -TMath::Log(TMath::Tan(theta / 2));
}

      

//____________________________________________________________________
void 
AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
			     TTree*  /* clusterTree */,
			     AliESDEvent* esd) const
{
  // nothing to be done
  // FIXME: The vertex may not be known when Reconstruct is executed,
  // so we may have to move some of that member function here. 
  AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
  // fESDObj->Print();

  if (esd) { 
    AliFMDDebug(2, ("Writing FMD data to ESD tree"));
    esd->SetFMDData(fESDObj);
  }

  if (!fDiagnostics || !esd) return;
  static bool first = true;
  // This is most likely NOT the event number you'd like to use. It
  // has nothing to do with the 'real' event number. 
  // - That's OK.  We just use it for the name of the directory -
  // nothing else.  Christian
  Int_t evno = esd->GetEventNumberInFile(); 
  AliFMDDebug(1, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
  TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
  first = false;
  f.cd(); 
  TDirectory* d = f.mkdir(Form("%03d", evno), 
			  Form("Diagnostics histograms for event # %d", evno));
  d->cd();
  if (fDiagStep1) fDiagStep1->Write();
  if (fDiagStep2) fDiagStep2->Write();
  if (fDiagStep3) fDiagStep3->Write();
  if (fDiagStep4) fDiagStep4->Write();
  if (fDiagAll)   fDiagAll->Write();
  d->Write();
  f.Write();
  f.Close();

  if (fDiagStep1) fDiagStep1->Reset();
  if (fDiagStep2) fDiagStep2->Reset();
  if (fDiagStep3) fDiagStep3->Reset();
  if (fDiagStep4) fDiagStep4->Reset();
  if (fDiagAll)   fDiagAll->Reset();
}

//____________________________________________________________________
void
AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree, 
			     AliESDEvent* esd) const
{
  TTree* dummy = 0;
  FillESD(dummy, clusterTree, esd);
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
37c55dc0	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	**************************************************************************/
0d0e6995	15	/* $Id$ */
c2fc1258	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	*/
4347b38f	21	//____________________________________________________________________
42403906	22	//
02a27b50	23	// This is a class that constructs AliFMDRecPoint objects from of Digits
6169f936	24	// This class reads either digits from a TClonesArray or raw data from
4347b38f	25	// a DDL file (or similar), and stores the read ADC counts in an
6169f936	26	// internal cache (fAdcs). The rec-points are made via the naiive
6169f936	27	// method.
4347b38f	28	//
37c55dc0	29	//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
4347b38f	30	// Latest changes by Christian Holm Christensen <cholm@nbi.dk>
	31	//
	32	//
	33	//____________________________________________________________________
37c55dc0	34
f95a63c4	35	// #include <AliLog.h> // ALILOG_H
6169f936	36	// #include <AliRun.h> // ALIRUN_H
f95a63c4	37	#include "AliFMDDebug.h"
1a1fdef7	38	#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
1a1fdef7	39	#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
56b1929b	40	#include "AliFMDDigit.h" // ALIFMDDIGIT_H
56b1929b	41	#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
56b1929b	42	#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
bf000c32	43	#include "AliFMDRecPoint.h" // ALIFMDMULTNAIIVE_H
af885e0f	44	#include "AliESDEvent.h" // ALIESDEVENT_H
aad72f45	45	#include "AliESDVertex.h" // ALIESDVERTEX_H
8f6ee336	46	#include <AliESDFMD.h> // ALIESDFMD_H
aad72f45	47	#include <TMath.h>
9684be2f	48	#include <TH1.h>
	49	#include <TH2.h>
	50	#include <TFile.h>
02a27b50	51	class AliRawReader;
4347b38f	52
4347b38f	53	//____________________________________________________________________
925e6570	54	ClassImp(AliFMDReconstructor)
1a1fdef7	55	#if 0
	56	; // This is here to keep Emacs for indenting the next line
	57	#endif
37c55dc0	58
4347b38f	59	//____________________________________________________________________
	60	AliFMDReconstructor::AliFMDReconstructor()
	61	: AliReconstructor(),
e9fd1e20	62	fMult(0x0),
	63	fNMult(0),
	64	fTreeR(0x0),
	65	fCurrentVertex(0),
	66	fESDObj(0x0),
9684be2f	67	fNoiseFactor(0),
	68	fAngleCorrect(kTRUE),
	69	fVertexType(kNoVertex),
9684be2f	70	fESD(0x0),
	71	fDiagnostics(kFALSE),
	72	fDiagStep1(0),
	73	fDiagStep2(0),
	74	fDiagStep3(0),
	75	fDiagStep4(0),
	76	fDiagAll(0)
4347b38f	77	{
8f6ee336	78	// Make a new FMD reconstructor object - default CTOR.
a9579262	79	SetNoiseFactor();
a9579262	80	SetAngleCorrect();
4347b38f	81	}
4347b38f	82
42403906	83
	84	//____________________________________________________________________
	85	AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
8f6ee336	86	: AliReconstructor(),
8f6ee336	87	fMult(other.fMult),
e9fd1e20	88	fNMult(other.fNMult),
	89	fTreeR(other.fTreeR),
	90	fCurrentVertex(other.fCurrentVertex),
	91	fESDObj(other.fESDObj),
a9579262	92	fNoiseFactor(other.fNoiseFactor),
9684be2f	93	fAngleCorrect(other.fAngleCorrect),
9684be2f	94	fVertexType(other.fVertexType),
9684be2f	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)
42403906	102	{
56b1929b	103	// Copy constructor
42403906	104	}
	105
	106
	107	//____________________________________________________________________
	108	AliFMDReconstructor&
	109	AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
	110	{
56b1929b	111	// Assignment operator
a9579262	112	fMult = other.fMult;
	113	fNMult = other.fNMult;
	114	fTreeR = other.fTreeR;
e9fd1e20	115	fCurrentVertex = other.fCurrentVertex;
a9579262	116	fESDObj = other.fESDObj;
a9579262	117	fNoiseFactor = other.fNoiseFactor;
a9579262	118	fAngleCorrect = other.fAngleCorrect;
9684be2f	119	fVertexType = other.fVertexType;
9684be2f	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;
42403906	127	return *this;
42403906	128	}
56b1929b	129
	130	//____________________________________________________________________
	131	AliFMDReconstructor::~AliFMDReconstructor()
	132	{
	133	// Destructor
8f6ee336	134	if (fMult) fMult->Delete();
	135	if (fMult) delete fMult;
	136	if (fESDObj) delete fESDObj;
56b1929b	137	}
4347b38f	138
	139	//____________________________________________________________________
	140	void
d76c31f4	141	AliFMDReconstructor::Init()
1a1fdef7	142	{
1a1fdef7	143	// Initialize the reconstructor
f011bd38	144
bf000c32	145	// Initialize the geometry
9b48326f	146	AliFMDGeometry* geom = AliFMDGeometry::Instance();
	147	geom->Init();
	148	geom->InitTransformations();
	149
	150	// Initialize the parameters
	151	AliFMDParameters* param = AliFMDParameters::Instance();
	152	param->Init();
bf000c32	153
8f6ee336	154	// Current vertex position
1a1fdef7	155	fCurrentVertex = 0;
8f6ee336	156	// Create array of reconstructed strip multiplicities
bf000c32	157	fMult = new TClonesArray("AliFMDRecPoint", 51200);
8f6ee336	158	// Create ESD output object
	159	fESDObj = new AliESDFMD;
	160
9684be2f	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");
4347b38f	189	}
dc8af42e	190
4347b38f	191	//____________________________________________________________________
4347b38f	192	void
1a1fdef7	193	AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
	194	TTree* digitsTree) const
	195	{
	196	// Convert Raw digits to AliFMDDigit's in a tree
f95a63c4	197	AliFMDDebug(2, ("Reading raw data into digits tree"));
1a1fdef7	198	AliFMDRawReader rawRead(reader, digitsTree);
	199	// rawRead.SetSampleRate(fFMD->GetSampleRate());
	200	rawRead.Exec();
4347b38f	201	}
4347b38f	202
4347b38f	203	//____________________________________________________________________
4347b38f	204	void
9684be2f	205	AliFMDReconstructor::GetVertex() const
4347b38f	206	{
97e94238	207	// Return the vertex to use.
	208	// This is obtained from the ESD object.
	209	// If not found, a warning is issued.
9684be2f	210	fVertexType = kNoVertex;
9684be2f	211	fCurrentVertex = 0;
a3537838	212	if (fESD) {
8f6ee336	213	const AliESDVertex* vertex = fESD->GetVertex();
8f6ee336	214	if (vertex) {
f95a63c4	215	AliFMDDebug(2, ("Got vertex from ESD: %f", vertex->GetZv()));
8f6ee336	216	fCurrentVertex = vertex->GetZv();
9684be2f	217	fVertexType = kESDVertex;
9684be2f	218	return;
8f6ee336	219	}
a3537838	220	}
9684be2f	221	AliWarning("Didn't get any vertex from ESD or generator");
	222	}
	223
	224
ddaa8027	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
9684be2f	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.
f95a63c4	254	AliFMDDebug(2, ("Reconstructing from digits in a tree"));
9684be2f	255	GetVertex();
9684be2f	256
1a1fdef7	257	TBranch *digitBranch = digitsTree->GetBranch("FMD");
4347b38f	258	if (!digitBranch) {
1a1fdef7	259	Error("Exec", "No digit branch for the FMD found");
1a1fdef7	260	return;
4347b38f	261	}
0abe7182	262	TClonesArray* digits = new TClonesArray("AliFMDDigit");
1a1fdef7	263	digitBranch->SetAddress(&digits);
4347b38f	264
8f6ee336	265	if (fMult) fMult->Clear();
	266	if (fESDObj) fESDObj->Clear();
	267
	268	fNMult = 0;
	269	fTreeR = clusterTree;
	270	fTreeR->Branch("FMD", &fMult);
	271
f95a63c4	272	AliFMDDebug(5, ("Getting entry 0 from digit branch"));
1a1fdef7	273	digitBranch->GetEntry(0);
1a1fdef7	274
f95a63c4	275	AliFMDDebug(5, ("Processing digits"));
e802be3e	276	ProcessDigits(digits);
e802be3e	277
8f6ee336	278	Int_t written = clusterTree->Fill();
f95a63c4	279	AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
0abe7182	280	digits->Delete();
0abe7182	281	delete digits;
4347b38f	282	}
1a1fdef7	283
8f6ee336	284
4347b38f	285	//____________________________________________________________________
e802be3e	286	void
e802be3e	287	AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f	288	{
69b696b9	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.
e802be3e	292	Int_t nDigits = digits->GetEntries();
f95a63c4	293	AliFMDDebug(1, ("Got %d digits", nDigits));
a9579262	294	fESDObj->SetNoiseFactor(fNoiseFactor);
a9579262	295	fESDObj->SetAngleCorrected(fAngleCorrect);
e802be3e	296	for (Int_t i = 0; i < nDigits; i++) {
e802be3e	297	AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
8f6ee336	298	AliFMDParameters* param = AliFMDParameters::Instance();
	299	// Check that the strip is not marked as dead
	300	if (param->IsDead(digit->Detector(), digit->Ring(),
15b17c89	301	digit->Sector(), digit->Strip())) {
f95a63c4	302	AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", digit->Detector(),
15b17c89	303	digit->Ring(), digit->Sector(), digit->Strip()));
	304	continue;
	305	}
8f6ee336	306
	307	// digit->Print();
	308	// Get eta and phi
	309	Float_t eta, phi;
	310	PhysicalCoordinates(digit, eta, phi);
4347b38f	311
8f6ee336	312	// Substract pedestal.
e802be3e	313	UShort_t counts = SubtractPedestal(digit);
4347b38f	314
8f6ee336	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
f95a63c4	321	AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
8f6ee336	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.
69893a66	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
8f6ee336	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);
9684be2f	343
9684be2f	344	if (fDiagAll) fDiagAll->Fill(digit->Counts(), mult);
4347b38f	345	}
4347b38f	346	}
8f6ee336	347
1a1fdef7	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
8f6ee336	357	Int_t counts = 0;
9684be2f	358	Int_t adc = 0;
8f6ee336	359	AliFMDParameters* param = AliFMDParameters::Instance();
a9579262	360	Float_t ped = param->GetPedestal(digit->Detector(),
8f6ee336	361	digit->Ring(),
	362	digit->Sector(),
	363	digit->Strip());
a9579262	364	Float_t noise = param->GetPedestalWidth(digit->Detector(),
	365	digit->Ring(),
	366	digit->Sector(),
	367	digit->Strip());
f95a63c4	368	AliFMDDebug(15, ("Subtracting pedestal %f from signal %d",
a9579262	369	ped, digit->Counts()));
9684be2f	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);
a9579262	374	if (counts < noise * fNoiseFactor) counts = 0;
f95a63c4	375	if (counts > 0) AliFMDDebug(15, ("Got a hit strip"));
9684be2f	376	if (fDiagStep1) fDiagStep1->Fill(adc, counts);
8f6ee336	377
1a1fdef7	378	return UShort_t(counts);
	379	}
	380
8f6ee336	381	//____________________________________________________________________
	382	Float_t
	383	AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit,
a9579262	384	Float_t eta,
8f6ee336	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.
a9579262	403	if (count <= 0) return 0;
8f6ee336	404	AliFMDParameters* param = AliFMDParameters::Instance();
	405	Float_t gain = param->GetPulseGain(digit->Detector(),
	406	digit->Ring(),
	407	digit->Sector(),
	408	digit->Strip());
f95a63c4	409	AliFMDDebug(15, ("Converting counts %d to energy via factor %f",
8f6ee336	410	count, gain));
a9579262	411
a9579262	412	Double_t edep = count * gain;
9684be2f	413	if (fDiagStep2) fDiagStep2->Fill(count, edep);
a9579262	414	if (fAngleCorrect) {
9684be2f	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;
f95a63c4	418	AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
9684be2f	419	edep, corr, cedep, eta, theta));
	420	if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
	421	edep = cedep;
a9579262	422	}
8f6ee336	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)
f95a63c4	449	AliFMDDebug(15, ("Translating energy %f to multiplicity via "
8f6ee336	450	"divider %f->%f", edep, edepMIP, mult));
9684be2f	451	if (fDiagStep4) fDiagStep4->Fill(edep, mult);
8f6ee336	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.
9b48326f	464	AliFMDGeometry* geom = AliFMDGeometry::Instance();
bf000c32	465	Double_t x, y, z, r, theta;
9b48326f	466	geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(),
bf000c32	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));
8f6ee336	474	}
	475
	476
	477
4347b38f	478	//____________________________________________________________________
4347b38f	479	void
1a1fdef7	480	AliFMDReconstructor::FillESD(TTree* /* digitsTree */,
1a1fdef7	481	TTree* /* clusterTree */,
af885e0f	482	AliESDEvent* esd) const
121a60bd	483	{
42403906	484	// nothing to be done
69b696b9	485	// FIXME: The vertex may not be known when Reconstruct is executed,
69b696b9	486	// so we may have to move some of that member function here.
f95a63c4	487	AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
8f6ee336	488	// fESDObj->Print();
	489
	490	if (esd) {
f95a63c4	491	AliFMDDebug(2, ("Writing FMD data to ESD tree"));
8f6ee336	492	esd->SetFMDData(fESDObj);
a3537838	493	}
9684be2f	494
	495	if (!fDiagnostics \|\| !esd) return;
	496	static bool first = true;
69893a66	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));
9684be2f	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();
121a60bd	523	}
121a60bd	524
ddaa8027	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
42403906	534	//____________________________________________________________________
	535	//
	536	// EOF
	537	//