[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$ */

//____________________________________________________________________
//
// This is a class that constructs ReconstParticles (reconstructed
// particles) out 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). 
//
// From the cached values it then calculates the number of particles
// that hit a region of the FMDs, as specified by the user. 
//
// The reconstruction can be done in two ways: Either via counting the
// number of empty strips (Poisson method), or by converting the ADC
// signal to an energy deposition, and then dividing by the typical
// energy loss of a particle.
// 
// Currently, this class only reads the digits from a TClonesArray,
// and the Poission method for reconstruction. 
//
// 
//-- 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 <AliRunLoader.h>                  // ALIRUNLOADER_H
#include <AliLoader.h>                     // ALILOADER_H
#include <AliHeader.h>                     // ALIHEADER_H
#include <AliRawReader.h>                  // ALIRAWREADER_H
#include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
#include "AliFMD.h"                        // ALIFMD_H
#include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
#include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
#include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
#include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
#include "AliFMDMultAlgorithm.h" 	   // ALIFMDMULTALGORITHM_H
#include "AliFMDMultPoisson.h"		   // ALIFMDMULTPOISSON_H
#include "AliFMDMultNaiive.h"		   // ALIFMDMULTNAIIVE_H

//____________________________________________________________________
ClassImp(AliFMDReconstructor);

//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor() 
  : AliReconstructor(),
    fPedestal(0), 
    fPedestalWidth(0),
    fPedestalFactor(0)
{
  // Make a new FMD reconstructor object - default CTOR.
  SetPedestal();

  fFMDLoader = 0;
  fRunLoader = 0;
  fFMD       = 0;
  fAlgorithms.Add(new AliFMDMultNaiive);
  fAlgorithms.Add(new AliFMDMultPoisson);
}
  

//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
  : AliReconstructor(),
    fPedestal(0), 
    fPedestalWidth(0),
    fPedestalFactor(0)
{
  // Copy constructor 
  SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);

  fFMDLoader = other.fFMDLoader;
  fRunLoader = other.fRunLoader;
  fFMD       = other.fFMD;
  
  fAlgorithms.Delete();
  TIter next(&(other.fAlgorithms));
  AliFMDMultAlgorithm* algorithm = 0;
  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
    fAlgorithms.Add(algorithm);
  fAlgorithms.SetOwner(kFALSE);
}
  

//____________________________________________________________________
AliFMDReconstructor&
AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
{
  // Assignment operator
  SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);

  fFMDLoader = other.fFMDLoader;
  fRunLoader = other.fRunLoader;
  fFMD       = other.fFMD;

  fAlgorithms.Delete();
  TIter next(&(other.fAlgorithms));
  AliFMDMultAlgorithm* algorithm = 0;
  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
    fAlgorithms.Add(algorithm);
  fAlgorithms.SetOwner(kFALSE);

  return *this;
}

//____________________________________________________________________
AliFMDReconstructor::~AliFMDReconstructor() 
{
  // Destructor 
  fAlgorithms.Delete();
}
  
//____________________________________________________________________
void 
AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor) 
{
  // Set the pedestal, and pedestal width 
  fPedestal       = mean;
  fPedestalWidth  = width;
  fPedestalFactor = factor;
}

//____________________________________________________________________
void 
AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader, 
				 AliRawReader* rawReader) const
{ 
  // Collects all digits in the same active volume into number of
  // particles
  //
  // Reconstruct number of particles in given group of pads for given
  // FMDvolume determined by numberOfVolume,
  // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
  // numberOgMaxRing 
  //
  // The reconstruction method is choosen based on the number of empty
  // strips. 
  if (!runLoader) {
    Error("Exec","Run Loader loader is NULL - Session not opened");
    return;
  }
  fRunLoader = runLoader;
  fFMDLoader = runLoader->GetLoader("FMDLoader");
  if (!fFMDLoader) 
    Fatal("AliFMDReconstructor","Can not find FMD (loader) "
	  "in specified event");

  // Get the AliRun object
  if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
  
  // Get the AliFMD object
  fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
  if (!fFMD) {
    AliError("Can not get FMD from gAlice");
    return;
  }
  fFMDLoader->LoadRecPoints("RECREATE");

  if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();

  TIter next(&fAlgorithms);
  AliFMDMultAlgorithm* algorithm = 0;
  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
    algorithm->PreRun(fFMD);

  if (rawReader) {
    Int_t event = 0;
    while (rawReader->NextEvent()) {
      ProcessEvent(event, rawReader);
      event++;
    }
  }
  else {
    Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries()); 
    for(Int_t event = 0; event < nEvents; event++) 
      ProcessEvent(event, 0);
  }

  next.Reset();
  algorithm = 0;
  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
    algorithm->PostRun();

  fFMDLoader->UnloadRecPoints();
  fFMDLoader = 0;
  fRunLoader = 0;
  fFMD       = 0;
}

//____________________________________________________________________
void 
AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
{ 
  // Collects all digits in the same active volume into number of
  // particles
  //
  // Reconstruct number of particles in given group of pads for given
  // FMDvolume determined by numberOfVolume,
  // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
  // numberOgMaxRing 
  //
  // The reconstruction method is choosen based on the number of empty
  // strips. 
  Reconstruct(runLoader, 0);
}


//____________________________________________________________________
void 
AliFMDReconstructor::ProcessEvent(Int_t event, 
				  AliRawReader* reader) const
{
  // Process one event read from either a clones array or from a a raw
  // data reader. 
  fRunLoader->GetEvent(event) ;
  //event z-vertex for correction eta-rad dependence      
  AliHeader *header            = fRunLoader->GetHeader();
  if (!header) Warning("ProcessEvent", "no AliHeader found!");
  AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);

  // Get the Z--coordinate from the event header 
  TArrayF o(3); 
  if (genHeader) genHeader->PrimaryVertex(o);
  else           Warning("ProcessEvent", "No GenEventHeader Found");
  fCurrentVertex = o.At(2);

  // If the recontruction tree isn't loaded, load it
  if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
  
  // Load or recreate the digits 
  if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
    if (!reader) {
      Error("Exec","Error occured while loading digits. Exiting.");
      return;
    }
  }
  // Get the digits tree 
  TTree* digitTree = fFMDLoader->TreeD();
  if (!digitTree) { 
    if (!reader) {
      Error("Exec","Can not get Tree with Digits. "
	    "Nothing to reconstruct - Exiting");
      return;
    }
    fFMDLoader->MakeTree("D");
    digitTree = fFMDLoader->TreeD();
    
  }
  // Get the FMD branch holding the digits. 
  TBranch *digitBranch = digitTree->GetBranch("FMD");
  TClonesArray* digits = fFMD->Digits();
  if (!digitBranch) {
    if (!reader) {
      Error("Exec", "No digit branch for the FMD found");
      return;
    }
    fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
  }
  if (!reader) digitBranch->SetAddress(&digits);

  if  (reader) {
    AliFMDRawReader rawRead(fFMD, reader);
    // rawRead->SetSampleRate(fSampleRate);
    rawRead.Exec();
  }
  else {
    // Read the ADC values from a clones array. 
    AliDebug(10, "Reading ADCs from Digits array");
    // read Digits, and reconstruct the particles
    if (!fFMDLoader->TreeD()->GetEvent(0)) return;
  }
  
  TIter next(&fAlgorithms);
  AliFMDMultAlgorithm* algorithm = 0;
  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
    algorithm->PreEvent(fFMDLoader->TreeR(), fCurrentVertex);

  ProcessDigits(digits);

  next.Reset();
  algorithm = 0;
  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
    algorithm->PostEvent();
  
  if (reader) {
    digitTree->Fill();
    fFMDLoader->WriteDigits("OVERWRITE");
  }
  fFMDLoader->UnloadDigits();
  fFMDLoader->TreeR()->Reset();
  fFMDLoader->TreeR()->Fill(); 
  fFMDLoader->WriteRecPoints("OVERWRITE");
}

//____________________________________________________________________
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;
  Float_t ped = fPedestal * fPedestalFactor * fPedestalWidth;
  if (digit->Count3() >= 0)      counts = digit->Count3();
  else if (digit->Count2() >= 0) counts = digit->Count2();
  else                           counts = digit->Count2();
  counts = TMath::Max(Int_t(counts - ped), 0);
  return  UShort_t(counts);
}

//____________________________________________________________________
void
AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
{
  Int_t nDigits = digits->GetEntries();
  for (Int_t i = 0; i < nDigits; i++) {
    AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
    AliFMDSubDetector* subDetector = 0;
    switch (digit->Detector()) {
    case 1: subDetector = fFMD->GetFMD1(); break;
    case 2: subDetector = fFMD->GetFMD2(); break;
    case 3: subDetector = fFMD->GetFMD3(); break;
    }
    if (!subDetector) { 
      Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
      continue;
    }
    
    AliFMDRing* ring  = 0;
    Float_t     ringZ = 0;
    switch(digit->Ring()) {
    case 'i':
    case 'I':  
      ring  = subDetector->GetInner(); 
      ringZ = subDetector->GetInnerZ();
      break;
    case 'o':
    case 'O':  
      ring  = subDetector->GetOuter(); 
      ringZ = subDetector->GetOuterZ();
      break;
    }
    if (!ring) {
      Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
	      digit->Ring());
      break;
    }
    
    Float_t  realZ    = fCurrentVertex + ringZ;
    Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
			 / ring->GetNStrips() * (digit->Strip() + .5) 
			 + ring->GetLowR());
    Float_t  theta    = TMath::ATan2(stripR, realZ);
    Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
			 * (digit->Sector() + .5));
    Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
    UShort_t counts   = SubtractPedestal(digit);
    
    TIter next(&fAlgorithms);
    AliFMDMultAlgorithm* algorithm = 0;
    while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
      algorithm->ProcessDigit(digit, eta, phi, counts);
  }
}
      
 
//____________________________________________________________________
void 
AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/, 
			     AliESD* /*esd*/) const
{
  // nothing to be done

}

//____________________________________________________________________
//
// 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	**************************************************************************/
	15
0d0e6995	16	/* $Id$ */
0d0e6995	17
4347b38f	18	//____________________________________________________________________
42403906	19	//
4347b38f	20	// This is a class that constructs ReconstParticles (reconstructed
	21	// particles) out of Digits
	22	//
	23	//
	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).
	27	//
	28	// From the cached values it then calculates the number of particles
	29	// that hit a region of the FMDs, as specified by the user.
	30	//
	31	// The reconstruction can be done in two ways: Either via counting the
	32	// number of empty strips (Poisson method), or by converting the ADC
	33	// signal to an energy deposition, and then dividing by the typical
	34	// energy loss of a particle.
	35	//
	36	// Currently, this class only reads the digits from a TClonesArray,
	37	// and the Poission method for reconstruction.
	38	//
37c55dc0	39	//
37c55dc0	40	//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
4347b38f	41	// Latest changes by Christian Holm Christensen <cholm@nbi.dk>
	42	//
	43	//
	44	//____________________________________________________________________
37c55dc0	45
56b1929b	46	#include <AliLog.h> // ALILOG_H
	47	#include <AliRun.h> // ALIRUN_H
	48	#include <AliRunLoader.h> // ALIRUNLOADER_H
	49	#include <AliLoader.h> // ALILOADER_H
	50	#include <AliHeader.h> // ALIHEADER_H
	51	#include <AliRawReader.h> // ALIRAWREADER_H
	52	#include <AliGenEventHeader.h> // ALIGENEVENTHEADER_H
	53	#include "AliFMD.h" // ALIFMD_H
	54	#include "AliFMDDigit.h" // ALIFMDDIGIT_H
	55	#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
	56	#include "AliFMDRawStream.h" // ALIFMDRAWSTREAM_H
	57	#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
	58	#include "AliFMDMultAlgorithm.h" // ALIFMDMULTALGORITHM_H
	59	#include "AliFMDMultPoisson.h" // ALIFMDMULTPOISSON_H
	60	#include "AliFMDMultNaiive.h" // ALIFMDMULTNAIIVE_H
4347b38f	61
	62	//____________________________________________________________________
	63	ClassImp(AliFMDReconstructor);
37c55dc0	64
4347b38f	65	//____________________________________________________________________
	66	AliFMDReconstructor::AliFMDReconstructor()
	67	: AliReconstructor(),
4347b38f	68	fPedestal(0),
e802be3e	69	fPedestalWidth(0),
e802be3e	70	fPedestalFactor(0)
4347b38f	71	{
42403906	72	// Make a new FMD reconstructor object - default CTOR.
4347b38f	73	SetPedestal();
4347b38f	74
4347b38f	75	fFMDLoader = 0;
	76	fRunLoader = 0;
	77	fFMD = 0;
56b1929b	78	fAlgorithms.Add(new AliFMDMultNaiive);
56b1929b	79	fAlgorithms.Add(new AliFMDMultPoisson);
4347b38f	80	}
4347b38f	81
42403906	82
	83	//____________________________________________________________________
	84	AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
	85	: AliReconstructor(),
42403906	86	fPedestal(0),
e802be3e	87	fPedestalWidth(0),
e802be3e	88	fPedestalFactor(0)
42403906	89	{
56b1929b	90	// Copy constructor
e802be3e	91	SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
42403906	92
42403906	93	fFMDLoader = other.fFMDLoader;
	94	fRunLoader = other.fRunLoader;
	95	fFMD = other.fFMD;
56b1929b	96
	97	fAlgorithms.Delete();
	98	TIter next(&(other.fAlgorithms));
	99	AliFMDMultAlgorithm* algorithm = 0;
	100	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
	101	fAlgorithms.Add(algorithm);
	102	fAlgorithms.SetOwner(kFALSE);
42403906	103	}
	104
	105
	106	//____________________________________________________________________
	107	AliFMDReconstructor&
	108	AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
	109	{
56b1929b	110	// Assignment operator
56b1929b	111	SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
42403906	112
42403906	113	fFMDLoader = other.fFMDLoader;
	114	fRunLoader = other.fRunLoader;
	115	fFMD = other.fFMD;
	116
56b1929b	117	fAlgorithms.Delete();
	118	TIter next(&(other.fAlgorithms));
	119	AliFMDMultAlgorithm* algorithm = 0;
	120	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
	121	fAlgorithms.Add(algorithm);
	122	fAlgorithms.SetOwner(kFALSE);
	123
42403906	124	return *this;
42403906	125	}
56b1929b	126
	127	//____________________________________________________________________
	128	AliFMDReconstructor::~AliFMDReconstructor()
	129	{
	130	// Destructor
	131	fAlgorithms.Delete();
	132	}
42403906	133
4347b38f	134	//____________________________________________________________________
4347b38f	135	void
e802be3e	136	AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor)
4347b38f	137	{
4347b38f	138	// Set the pedestal, and pedestal width
e802be3e	139	fPedestal = mean;
	140	fPedestalWidth = width;
	141	fPedestalFactor = factor;
4347b38f	142	}
	143
	144	//____________________________________________________________________
	145	void
	146	AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
	147	AliRawReader* rawReader) const
	148	{
	149	// Collects all digits in the same active volume into number of
	150	// particles
	151	//
	152	// Reconstruct number of particles in given group of pads for given
	153	// FMDvolume determined by numberOfVolume,
	154	// numberOfMinSector, numberOfMaxSector, numberOfMinRing,
	155	// numberOgMaxRing
	156	//
	157	// The reconstruction method is choosen based on the number of empty
	158	// strips.
4347b38f	159	if (!runLoader) {
	160	Error("Exec","Run Loader loader is NULL - Session not opened");
	161	return;
	162	}
	163	fRunLoader = runLoader;
	164	fFMDLoader = runLoader->GetLoader("FMDLoader");
	165	if (!fFMDLoader)
	166	Fatal("AliFMDReconstructor","Can not find FMD (loader) "
	167	"in specified event");
dc8af42e	168
4347b38f	169	// Get the AliRun object
	170	if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
	171
	172	// Get the AliFMD object
	173	fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
	174	if (!fFMD) {
	175	AliError("Can not get FMD from gAlice");
	176	return;
	177	}
	178	fFMDLoader->LoadRecPoints("RECREATE");
dc8af42e	179
4347b38f	180	if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
88cb7938	181
56b1929b	182	TIter next(&fAlgorithms);
	183	AliFMDMultAlgorithm* algorithm = 0;
	184	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
	185	algorithm->PreRun(fFMD);
	186
4347b38f	187	if (rawReader) {
	188	Int_t event = 0;
	189	while (rawReader->NextEvent()) {
e802be3e	190	ProcessEvent(event, rawReader);
4347b38f	191	event++;
	192	}
	193	}
	194	else {
	195	Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
	196	for(Int_t event = 0; event < nEvents; event++)
e802be3e	197	ProcessEvent(event, 0);
4347b38f	198	}
88cb7938	199
56b1929b	200	next.Reset();
	201	algorithm = 0;
	202	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
	203	algorithm->PostRun();
dc8af42e	204
4347b38f	205	fFMDLoader->UnloadRecPoints();
	206	fFMDLoader = 0;
	207	fRunLoader = 0;
	208	fFMD = 0;
	209	}
dc8af42e	210
4347b38f	211	//____________________________________________________________________
	212	void
	213	AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e	214	{
4347b38f	215	// Collects all digits in the same active volume into number of
	216	// particles
	217	//
	218	// Reconstruct number of particles in given group of pads for given
	219	// FMDvolume determined by numberOfVolume,
	220	// numberOfMinSector, numberOfMaxSector, numberOfMinRing,
	221	// numberOgMaxRing
	222	//
	223	// The reconstruction method is choosen based on the number of empty
	224	// strips.
	225	Reconstruct(runLoader, 0);
	226	}
	227
	228
	229	//____________________________________________________________________
	230	void
	231	AliFMDReconstructor::ProcessEvent(Int_t event,
e802be3e	232	AliRawReader* reader) const
4347b38f	233	{
42403906	234	// Process one event read from either a clones array or from a a raw
42403906	235	// data reader.
4347b38f	236	fRunLoader->GetEvent(event) ;
	237	//event z-vertex for correction eta-rad dependence
	238	AliHeader *header = fRunLoader->GetHeader();
e802be3e	239	if (!header) Warning("ProcessEvent", "no AliHeader found!");
4347b38f	240	AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
3d44ce66	241
4347b38f	242	// Get the Z--coordinate from the event header
	243	TArrayF o(3);
	244	if (genHeader) genHeader->PrimaryVertex(o);
e802be3e	245	else Warning("ProcessEvent", "No GenEventHeader Found");
e802be3e	246	fCurrentVertex = o.At(2);
cb1df35e	247
4347b38f	248	// If the recontruction tree isn't loaded, load it
4347b38f	249	if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
88cb7938	250
4347b38f	251	// Load or recreate the digits
	252	if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
	253	if (!reader) {
	254	Error("Exec","Error occured while loading digits. Exiting.");
	255	return;
	256	}
4347b38f	257	}
	258	// Get the digits tree
	259	TTree* digitTree = fFMDLoader->TreeD();
	260	if (!digitTree) {
	261	if (!reader) {
	262	Error("Exec","Can not get Tree with Digits. "
	263	"Nothing to reconstruct - Exiting");
	264	return;
	265	}
	266	fFMDLoader->MakeTree("D");
	267	digitTree = fFMDLoader->TreeD();
	268
	269	}
	270	// Get the FMD branch holding the digits.
	271	TBranch *digitBranch = digitTree->GetBranch("FMD");
e802be3e	272	TClonesArray* digits = fFMD->Digits();
4347b38f	273	if (!digitBranch) {
	274	if (!reader) {
	275	Error("Exec", "No digit branch for the FMD found");
	276	return;
	277	}
	278	fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
	279	}
	280	if (!reader) digitBranch->SetAddress(&digits);
	281
e802be3e	282	if (reader) {
	283	AliFMDRawReader rawRead(fFMD, reader);
	284	// rawRead->SetSampleRate(fSampleRate);
	285	rawRead.Exec();
	286	}
	287	else {
	288	// Read the ADC values from a clones array.
	289	AliDebug(10, "Reading ADCs from Digits array");
	290	// read Digits, and reconstruct the particles
	291	if (!fFMDLoader->TreeD()->GetEvent(0)) return;
	292	}
4347b38f	293
e802be3e	294	TIter next(&fAlgorithms);
56b1929b	295	AliFMDMultAlgorithm* algorithm = 0;
	296	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
	297	algorithm->PreEvent(fFMDLoader->TreeR(), fCurrentVertex);
4347b38f	298
e802be3e	299	ProcessDigits(digits);
	300
	301	next.Reset();
	302	algorithm = 0;
56b1929b	303	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
e802be3e	304	algorithm->PostEvent();
e802be3e	305
4347b38f	306	if (reader) {
	307	digitTree->Fill();
	308	fFMDLoader->WriteDigits("OVERWRITE");
	309	}
	310	fFMDLoader->UnloadDigits();
	311	fFMDLoader->TreeR()->Reset();
	312	fFMDLoader->TreeR()->Fill();
	313	fFMDLoader->WriteRecPoints("OVERWRITE");
	314	}
	315
	316	//____________________________________________________________________
e802be3e	317	UShort_t
e802be3e	318	AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
4347b38f	319	{
e802be3e	320	// Member function to subtract the pedestal from a digit
	321	// This implementation does nothing, but a derived class could over
	322	// load this to subtract a pedestal that was given in a database or
	323	// something like that.
4347b38f	324
e802be3e	325	Int_t counts = 0;
	326	Float_t ped = fPedestal * fPedestalFactor * fPedestalWidth;
	327	if (digit->Count3() >= 0) counts = digit->Count3();
	328	else if (digit->Count2() >= 0) counts = digit->Count2();
	329	else counts = digit->Count2();
	330	counts = TMath::Max(Int_t(counts - ped), 0);
	331	return UShort_t(counts);
4347b38f	332	}
	333
	334	//____________________________________________________________________
e802be3e	335	void
e802be3e	336	AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f	337	{
e802be3e	338	Int_t nDigits = digits->GetEntries();
	339	for (Int_t i = 0; i < nDigits; i++) {
	340	AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
	341	AliFMDSubDetector* subDetector = 0;
	342	switch (digit->Detector()) {
	343	case 1: subDetector = fFMD->GetFMD1(); break;
	344	case 2: subDetector = fFMD->GetFMD2(); break;
	345	case 3: subDetector = fFMD->GetFMD3(); break;
	346	}
	347	if (!subDetector) {
	348	Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
	349	continue;
	350	}
4347b38f	351
e802be3e	352	AliFMDRing* ring = 0;
	353	Float_t ringZ = 0;
	354	switch(digit->Ring()) {
	355	case 'i':
	356	case 'I':
	357	ring = subDetector->GetInner();
	358	ringZ = subDetector->GetInnerZ();
	359	break;
	360	case 'o':
	361	case 'O':
	362	ring = subDetector->GetOuter();
	363	ringZ = subDetector->GetOuterZ();
	364	break;
	365	}
	366	if (!ring) {
	367	Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
	368	digit->Ring());
	369	break;
4347b38f	370	}
4347b38f	371
e802be3e	372	Float_t realZ = fCurrentVertex + ringZ;
56b1929b	373	Float_t stripR = ((ring->GetHighR() - ring->GetLowR())
	374	/ ring->GetNStrips() * (digit->Strip() + .5)
	375	+ ring->GetLowR());
e802be3e	376	Float_t theta = TMath::ATan2(stripR, realZ);
	377	Float_t phi = (2 * TMath::Pi() / ring->GetNSectors()
	378	* (digit->Sector() + .5));
	379	Float_t eta = -TMath::Log(TMath::Tan(theta / 2));
	380	UShort_t counts = SubtractPedestal(digit);
4347b38f	381
e802be3e	382	TIter next(&fAlgorithms);
56b1929b	383	AliFMDMultAlgorithm* algorithm = 0;
56b1929b	384	while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
e802be3e	385	algorithm->ProcessDigit(digit, eta, phi, counts);
4347b38f	386	}
4347b38f	387	}
4347b38f	388
4347b38f	389
	390	//____________________________________________________________________
	391	void
	392	AliFMDReconstructor::FillESD(AliRunLoader* /fRunLoader/,
	393	AliESD* /esd/) const
121a60bd	394	{
42403906	395	// nothing to be done
121a60bd	396
	397	}
	398
42403906	399	//____________________________________________________________________
	400	//
	401	// EOF
	402	//