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

//____________________________________________________________________
// 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 "AliFMD.h"
#include "AliFMDDigit.h"
#include "AliFMDParticles.h"
#include "AliFMDReconstructor.h"
#include "AliAltroBuffer.h"
#include "AliLog.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliFMDRawStream.h"
#include "AliRawReader.h"

//____________________________________________________________________
ClassImp(AliFMDReconstructor);

//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor() 
  : AliReconstructor(),
    fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
    fDeltaEta(0), 
    fDeltaPhi(0), 
    fThreshold(0),
    fPedestal(0), 
    fPedestalWidth(0)
{
  SetDeltaEta();
  SetDeltaPhi();
  SetThreshold();
  SetPedestal();

  fParticles = new TClonesArray("AliFMDParticles", 1000);
  fFMDLoader = 0;
  fRunLoader = 0;
  fFMD       = 0;
}
  
//____________________________________________________________________
void 
AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width) 
{
  // Set the pedestal, and pedestal width 
  fPedestal      = mean;
  fPedestalWidth = width;
}

//____________________________________________________________________
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. 
  fParticles->Clear();
  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();

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


  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, 
				  TClonesArray* digits) const
{
  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);
  Float_t zVertex = o.At(2);

  // If the recontruction tree isn't loaded, load it
  if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
  
  //Make branches to hold the reconstructed particles 
  const Int_t kBufferSize = 16000;
  fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);

  // 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");
  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);

  fEmptyStrips = 0;
  fTotalStrips = 0;
  Bool_t ok = kFALSE;
  if      (reader) ok = ReadAdcs(reader);
  else if (digits) ok = ReadAdcs(digits);
  if (!ok) return;
  
  ReconstructFromCache(zVertex);

  if (reader) {
    digitTree->Fill();
    fFMDLoader->WriteDigits("OVERWRITE");
  }
  fFMDLoader->UnloadDigits();
  fFMDLoader->TreeR()->Reset();
  fFMDLoader->TreeR()->Fill(); 
  fFMDLoader->WriteRecPoints("OVERWRITE");
}

//____________________________________________________________________
Bool_t
AliFMDReconstructor::ReadAdcs(TClonesArray* digits) const
{
  AliDebug(10, "Reading ADCs from Digits array");
  // read Digits, and reconstruct the particles
  if (!fFMDLoader->TreeD()->GetEvent(0)) return kFALSE;

  // Reads the digits from the array, and fills up the cache (fAdcs) 
  fAdcs.Clear();
  Int_t nDigits = digits->GetEntries();
  for (Int_t digit = 0; digit < nDigits; digit++) {
    AliFMDDigit* fmdDigit = 
      static_cast<AliFMDDigit*>(digits->UncheckedAt(digit));    
    
    ProcessDigit(fmdDigit);
  } //digit loop
  return kTRUE;
}

//____________________________________________________________________
Bool_t
AliFMDReconstructor::ReadAdcs(AliRawReader* reader) const
{
  AliDebug(10, "Reading ADCs from RawReader");
  // Reads the digits from a RAW data 
  fAdcs.Clear();
  // reader->Reset();

  if (!reader->ReadHeader()) {
    Error("ReadAdcs", "Couldn't read header");
    return kFALSE;
  }

  // Use AliAltroRawStream to read the ALTRO format.  No need to
  // reinvent the wheel :-) 
  AliFMDRawStream input(reader);
  // Select FMD DDL's 
  reader->Select(AliFMD::kBaseDDL >> 8);
  
  Int_t    oldDDL      = -1;
  Int_t    count       = 0;
  UShort_t detector    = 1; // Must be one here
  UShort_t oldDetector = 0;
  // Loop over data in file 
  Bool_t   next       = kTRUE;

  // local Cache 
  TArrayI counts(10);
  counts.Reset(-1);
  Int_t offset = 0;
  
  while (next) {
    next = input.Next();


    count++; 
    Int_t ddl = reader->GetDDLID();
    if (ddl != oldDDL 
	|| input.IsNewStrip()
	|| !next) {
      // Make a new digit, if we have some data! 
      if (counts[0] >= 0) {
	// Got a new strip. 
	AliDebug(10, Form("Add a new strip: FMD%d%c[%2d,%3d] "
			  "(current: FMD%d%c[%2d,%3d])", 
			  oldDetector, input.PrevRing(), 
			  input.PrevSector() , input.PrevStrip(),
			  detector , input.Ring(), input.Sector(), 
			  input.Strip()));
	fFMD->AddDigit(oldDetector, 
		       input.PrevRing(), 
		       input.PrevSector(), 
		       input.PrevStrip(), 
		       counts[0], counts[1], counts[2]);
	AliFMDDigit* digit = 
	  static_cast<AliFMDDigit*>(fFMD->Digits()->
				    UncheckedAt(fFMD->GetNdigits()-1));
	ProcessDigit(digit);
      }
	
      if (!next) { 
	AliDebug(10, Form("Read %d channels for FMD%d", 
			  count + 1, detector));
	break;
      }
    
    
      // If we got a new DDL, it means we have a new detector. 
      if (ddl != oldDDL) {
	if (detector != 0) 
	  AliDebug(10, Form("Read %d channels for FMD%d", 
			    count + 1, detector));
	// Reset counts, and update the DDL cache 
	count       = 0;
	oldDDL      = ddl;
	// Check that we're processing a FMD detector 
	Int_t detId = reader->GetDetectorID();
	if (detId != (AliFMD::kBaseDDL >> 8)) {
	  Error("ReadAdcs", "Detector ID %d != %d",
		detId, (AliFMD::kBaseDDL >> 8));
	  break;
	}
	// Figure out what detector we're deling with 
	oldDetector = detector;
	switch (ddl) {
	case 0: detector = 1; break;
	case 1: detector = 2; break;
	case 2: detector = 3; break;
	default:
	  Error("ReadAdcs", "Unknown DDL 0x%x for FMD", ddl);
	  return kFALSE;
	}
	AliDebug(10, Form("Reading ADCs for 0x%x  - That is FMD%d",
			  reader->GetEquipmentId(), detector));
      }
      counts.Reset(-1);
      offset = 0;
    }
    
    counts[offset] = input.Count();
    offset++;
    
    AliDebug(10, Form("ADC of FMD%d%c[%2d,%3d] += %d",
		      detector, input.Ring(), input.Sector(), 
		      input.Strip(), input.Count()));
    oldDetector = detector;
  }
  return kTRUE;
}

//____________________________________________________________________
void
AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
{
  // Process a digit.  Derived classes can overload this member
  // function to do stuff to the digit.  However, it should write the
  // ADC count to the internal cache 
  //
  //    fAdcs(detector - 1, ring, sector, strip) = counts;
  //
  // In this implementation, we count the number of strips below
  // threshold.   This we do to later choose what kind of
  // reconstruction algorithm we'd like to use. 
  // 
  UShort_t detector = digit->Detector();
  Char_t   ring     = digit->Ring();
  UShort_t sector   = digit->Sector();
  UShort_t strip    = digit->Strip();    
  
  UShort_t counts  = SubtractPedestal(digit);
  
  fAdcs(detector - 1, ring, sector, strip) = counts;
  if (counts < fThreshold) fEmptyStrips++;
  fTotalStrips++;
}

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

//____________________________________________________________________
void
AliFMDReconstructor::ReconstructFromCache(Float_t zVertex) const
{
  // Based on the information in the cache, do the reconstruction. 
  Int_t nRecon = 0;
  // Loop over the detectors 
  for (Int_t i = 1; i <= 3; i++) {
    AliFMDSubDetector* sub = 0;
    switch (i) {
    case 1: sub = fFMD->GetFMD1(); break;
    case 2: sub = fFMD->GetFMD2(); break;
    case 3: sub = fFMD->GetFMD3(); break;
    }
    if (!sub) continue;
	
    // Loop over the rings in the detector
    for (Int_t j = 0; j < 2; j++) {
      Float_t     rZ = 0;
      AliFMDRing* r  = 0;
      switch (j) {
      case 0: r  = sub->GetInner(); rZ = sub->GetInnerZ(); break;
      case 1: r  = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
      }
      if (!r) continue;
      
      // Calculate low/high theta and eta 
      // FIXME: Is this right? 
      Float_t realZ    = zVertex + rZ;
      Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
      Float_t thetaIn  = TMath::ATan2(r->GetLowR(), realZ);
      Float_t etaOut   = - TMath::Log(TMath::Tan(thetaOut / 2));
      Float_t etaIn    = - TMath::Log(TMath::Tan(thetaIn / 2));
      if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
	Float_t tmp = etaIn;
	etaIn       = etaOut;
	etaOut      = tmp;
      }

      //-------------------------------------------------------------
      //
      // Here starts poisson method 
      //
      // Calculate eta step per strip, number of eta steps, number of
      // phi steps, and check the sign of the eta increment 
      Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
      Int_t   nEta     = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta); 
      Int_t   nPhi     = Int_t(360. / fDeltaPhi);
      Float_t sign     = TMath::Sign(Float_t(1.), etaIn);

      AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
			sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));

      // Loop over relevant phi values 
      for (Int_t p = 0; p < nPhi; p++) {
	Float_t  minPhi    = p * fDeltaPhi;
	Float_t  maxPhi    = minPhi + fDeltaPhi;
	UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
	UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
	
	AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
			  minPhi, maxPhi, minSector, maxSector));
	// Loop over relevant eta values 
	for (Int_t e = nEta; e >= 0; --e) {
	  Float_t  maxEta   = etaIn  - sign * e * fDeltaEta;
	  Float_t  minEta   = maxEta - sign * fDeltaEta;
	  if (sign > 0)  minEta = TMath::Max(minEta, etaOut);
	  else           minEta = TMath::Min(minEta, etaOut);
	  Float_t  theta1   = 2 * TMath::ATan(TMath::Exp(-minEta));
	  Float_t  theta2   = 2 * TMath::ATan(TMath::Exp(-maxEta));
	  Float_t  minR     = TMath::Abs(realZ * TMath::Tan(theta2));
	  Float_t  maxR     = TMath::Abs(realZ * TMath::Tan(theta1));
	  UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
	  UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);

	  AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
			    "    [radii %f, %f, thetas %f, %f, sign %d]", 
			    minEta, maxEta, minStrip, maxStrip, 
			    minR, maxR, theta1, theta2, sign));

	  // Count number of empty strips
	  Int_t   emptyStrips = 0;
	  for (Int_t sector = minSector; sector < maxSector; sector++) 
	    for (Int_t strip = minStrip; strip < maxStrip; strip++) 
	      if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip) 
		  < fThreshold) emptyStrips++;
	  
	  // The total number of strips 
	  Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
	  
	  // Log ratio of empty to total number of strips 
	  AliDebug(10, Form("Lambda= %d / %d = %f", 
			    emptyStrips, nTotal, 
			    Float_t(emptyStrips) / nTotal));
	  
	  Double_t lambda = (emptyStrips > 0 ? 
			     - TMath::Log(Double_t(emptyStrips) / nTotal) :
			     1);

	  // The reconstructed number of particles is then given by 
	  Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
	    
	  // Add a AliFMDParticles to the reconstruction tree. 
	  new((*fParticles)[nRecon])   
	    AliFMDParticles(sub->GetId(), r->GetId(),
			    minSector, maxSector, minStrip, maxStrip,
			    minEta, maxEta, minPhi, maxPhi,
			    reconstructed, AliFMDParticles::kPoission);
	  nRecon++;
	} // phi 
      } // eta
    } // ring 
  } // detector 
}

 
//____________________________________________________________________
void 
AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/, 
			     AliESD* /*esd*/) const
{
// nothing to be done

}
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
4347b38f	16	//____________________________________________________________________
	17	// This is a class that constructs ReconstParticles (reconstructed
	18	// particles) out of Digits
	19	//
	20	//
	21	// This class reads either digits from a TClonesArray or raw data from
	22	// a DDL file (or similar), and stores the read ADC counts in an
	23	// internal cache (fAdcs).
	24	//
	25	// From the cached values it then calculates the number of particles
	26	// that hit a region of the FMDs, as specified by the user.
	27	//
	28	// The reconstruction can be done in two ways: Either via counting the
	29	// number of empty strips (Poisson method), or by converting the ADC
	30	// signal to an energy deposition, and then dividing by the typical
	31	// energy loss of a particle.
	32	//
	33	// Currently, this class only reads the digits from a TClonesArray,
	34	// and the Poission method for reconstruction.
	35	//
37c55dc0	36	//
37c55dc0	37	//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
4347b38f	38	// Latest changes by Christian Holm Christensen <cholm@nbi.dk>
	39	//
	40	//
	41	//____________________________________________________________________
37c55dc0	42
8f1cfb0c	43	#include "AliFMD.h"
	44	#include "AliFMDDigit.h"
	45	#include "AliFMDParticles.h"
	46	#include "AliFMDReconstructor.h"
	47	#include "AliAltroBuffer.h"
	48	#include "AliLog.h"
	49	#include "AliRun.h"
	50	#include "AliRunLoader.h"
	51	#include "AliLoader.h"
	52	#include "AliHeader.h"
	53	#include "AliGenEventHeader.h"
	54	#include "AliFMDRawStream.h"
	55	#include "AliRawReader.h"
4347b38f	56
	57	//____________________________________________________________________
	58	ClassImp(AliFMDReconstructor);
37c55dc0	59
4347b38f	60	//____________________________________________________________________
	61	AliFMDReconstructor::AliFMDReconstructor()
	62	: AliReconstructor(),
	63	fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
	64	fDeltaEta(0),
	65	fDeltaPhi(0),
	66	fThreshold(0),
	67	fPedestal(0),
	68	fPedestalWidth(0)
	69	{
	70	SetDeltaEta();
	71	SetDeltaPhi();
	72	SetThreshold();
	73	SetPedestal();
	74
	75	fParticles = new TClonesArray("AliFMDParticles", 1000);
	76	fFMDLoader = 0;
	77	fRunLoader = 0;
	78	fFMD = 0;
	79	}
	80
	81	//____________________________________________________________________
	82	void
	83	AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width)
	84	{
	85	// Set the pedestal, and pedestal width
	86	fPedestal = mean;
	87	fPedestalWidth = width;
	88	}
	89
	90	//____________________________________________________________________
	91	void
	92	AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
	93	AliRawReader* rawReader) const
	94	{
	95	// Collects all digits in the same active volume into number of
	96	// particles
	97	//
	98	// Reconstruct number of particles in given group of pads for given
	99	// FMDvolume determined by numberOfVolume,
	100	// numberOfMinSector, numberOfMaxSector, numberOfMinRing,
	101	// numberOgMaxRing
	102	//
	103	// The reconstruction method is choosen based on the number of empty
	104	// strips.
	105	fParticles->Clear();
	106	if (!runLoader) {
	107	Error("Exec","Run Loader loader is NULL - Session not opened");
	108	return;
	109	}
	110	fRunLoader = runLoader;
	111	fFMDLoader = runLoader->GetLoader("FMDLoader");
	112	if (!fFMDLoader)
	113	Fatal("AliFMDReconstructor","Can not find FMD (loader) "
	114	"in specified event");
dc8af42e	115
4347b38f	116	// Get the AliRun object
	117	if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
	118
	119	// Get the AliFMD object
	120	fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
	121	if (!fFMD) {
	122	AliError("Can not get FMD from gAlice");
	123	return;
	124	}
	125	fFMDLoader->LoadRecPoints("RECREATE");
dc8af42e	126
4347b38f	127	if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
88cb7938	128
4347b38f	129	TClonesArray* digits = fFMD->Digits();
	130	if (rawReader) {
	131	Int_t event = 0;
	132	while (rawReader->NextEvent()) {
	133	ProcessEvent(event, rawReader, digits);
	134	event++;
	135	}
	136	}
	137	else {
	138	Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
	139	for(Int_t event = 0; event < nEvents; event++)
	140	ProcessEvent(event, 0, digits);
	141	}
88cb7938	142
dc8af42e	143
4347b38f	144	fFMDLoader->UnloadRecPoints();
	145	fFMDLoader = 0;
	146	fRunLoader = 0;
	147	fFMD = 0;
	148	}
dc8af42e	149
4347b38f	150	//____________________________________________________________________
	151	void
	152	AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e	153	{
4347b38f	154	// Collects all digits in the same active volume into number of
	155	// particles
	156	//
	157	// Reconstruct number of particles in given group of pads for given
	158	// FMDvolume determined by numberOfVolume,
	159	// numberOfMinSector, numberOfMaxSector, numberOfMinRing,
	160	// numberOgMaxRing
	161	//
	162	// The reconstruction method is choosen based on the number of empty
	163	// strips.
	164	Reconstruct(runLoader, 0);
	165	}
	166
	167
	168	//____________________________________________________________________
	169	void
	170	AliFMDReconstructor::ProcessEvent(Int_t event,
	171	AliRawReader* reader,
	172	TClonesArray* digits) const
	173	{
	174	fRunLoader->GetEvent(event) ;
	175	//event z-vertex for correction eta-rad dependence
	176	AliHeader *header = fRunLoader->GetHeader();
	177	if (!header)
	178	Warning("ProcessEvent", "no AliHeader found!");
	179	AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
3d44ce66	180
4347b38f	181	// Get the Z--coordinate from the event header
	182	TArrayF o(3);
	183	if (genHeader) genHeader->PrimaryVertex(o);
	184	Float_t zVertex = o.At(2);
cb1df35e	185
4347b38f	186	// If the recontruction tree isn't loaded, load it
4347b38f	187	if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
88cb7938	188
4347b38f	189	//Make branches to hold the reconstructed particles
	190	const Int_t kBufferSize = 16000;
	191	fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
	192
	193	// Load or recreate the digits
	194	if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
	195	if (!reader) {
	196	Error("Exec","Error occured while loading digits. Exiting.");
	197	return;
	198	}
46501dfb	199
4347b38f	200	}
	201	// Get the digits tree
	202	TTree* digitTree = fFMDLoader->TreeD();
	203	if (!digitTree) {
	204	if (!reader) {
	205	Error("Exec","Can not get Tree with Digits. "
	206	"Nothing to reconstruct - Exiting");
	207	return;
	208	}
	209	fFMDLoader->MakeTree("D");
	210	digitTree = fFMDLoader->TreeD();
	211
	212	}
	213	// Get the FMD branch holding the digits.
	214	TBranch *digitBranch = digitTree->GetBranch("FMD");
	215	if (!digitBranch) {
	216	if (!reader) {
	217	Error("Exec", "No digit branch for the FMD found");
	218	return;
	219	}
	220	fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
	221	}
	222	if (!reader) digitBranch->SetAddress(&digits);
	223
	224	fEmptyStrips = 0;
	225	fTotalStrips = 0;
	226	Bool_t ok = kFALSE;
	227	if (reader) ok = ReadAdcs(reader);
	228	else if (digits) ok = ReadAdcs(digits);
	229	if (!ok) return;
	230
	231	ReconstructFromCache(zVertex);
	232
	233	if (reader) {
	234	digitTree->Fill();
	235	fFMDLoader->WriteDigits("OVERWRITE");
	236	}
	237	fFMDLoader->UnloadDigits();
	238	fFMDLoader->TreeR()->Reset();
	239	fFMDLoader->TreeR()->Fill();
	240	fFMDLoader->WriteRecPoints("OVERWRITE");
	241	}
	242
	243	//____________________________________________________________________
	244	Bool_t
	245	AliFMDReconstructor::ReadAdcs(TClonesArray* digits) const
	246	{
	247	AliDebug(10, "Reading ADCs from Digits array");
	248	// read Digits, and reconstruct the particles
	249	if (!fFMDLoader->TreeD()->GetEvent(0)) return kFALSE;
	250
	251	// Reads the digits from the array, and fills up the cache (fAdcs)
	252	fAdcs.Clear();
	253	Int_t nDigits = digits->GetEntries();
	254	for (Int_t digit = 0; digit < nDigits; digit++) {
	255	AliFMDDigit* fmdDigit =
	256	static_cast<AliFMDDigit*>(digits->UncheckedAt(digit));
	257
	258	ProcessDigit(fmdDigit);
	259	} //digit loop
	260	return kTRUE;
	261	}
	262
	263	//____________________________________________________________________
264	Bool_t
265	AliFMDReconstructor::ReadAdcs(AliRawReader* reader) const
266	{
267	AliDebug(10, "Reading ADCs from RawReader");
268	// Reads the digits from a RAW data
269	fAdcs.Clear();
270	// reader->Reset();
271
272	if (!reader->ReadHeader()) {
273	Error("ReadAdcs", "Couldn't read header");
274	return kFALSE;
275	}
276
277	// Use AliAltroRawStream to read the ALTRO format. No need to
278	// reinvent the wheel :-)
279	AliFMDRawStream input(reader);
280	// Select FMD DDL's
281	reader->Select(AliFMD::kBaseDDL >> 8);
282
283	Int_t oldDDL = -1;
284	Int_t count = 0;
285	UShort_t detector = 1; // Must be one here
286	UShort_t oldDetector = 0;
287	// Loop over data in file
288	Bool_t next = kTRUE;
289
290	// local Cache
291	TArrayI counts(10);
292	counts.Reset(-1);
293	Int_t offset = 0;
294
295	while (next) {
296	next = input.Next();
297
298
299	count++;
300	Int_t ddl = reader->GetDDLID();
301	if (ddl != oldDDL
302	\|\| input.IsNewStrip()
303	\|\| !next) {
304	// Make a new digit, if we have some data!
305	if (counts[0] >= 0) {
306	// Got a new strip.
307	AliDebug(10, Form("Add a new strip: FMD%d%c[%2d,%3d] "
308	"(current: FMD%d%c[%2d,%3d])",
309	oldDetector, input.PrevRing(),
310	input.PrevSector() , input.PrevStrip(),
311	detector , input.Ring(), input.Sector(),
312	input.Strip()));
313	fFMD->AddDigit(oldDetector,
314	input.PrevRing(),
315	input.PrevSector(),
316	input.PrevStrip(),
317	counts[0], counts[1], counts[2]);
318	AliFMDDigit* digit =
319	static_cast<AliFMDDigit*>(fFMD->Digits()->
320	UncheckedAt(fFMD->GetNdigits()-1));
321	ProcessDigit(digit);
322	}
323
324	if (!next) {
325	AliDebug(10, Form("Read %d channels for FMD%d",
326	count + 1, detector));
327	break;
328	}
329
330
331	// If we got a new DDL, it means we have a new detector.
332	if (ddl != oldDDL) {
333	if (detector != 0)
334	AliDebug(10, Form("Read %d channels for FMD%d",
335	count + 1, detector));
336	// Reset counts, and update the DDL cache
337	count = 0;
338	oldDDL = ddl;
339	// Check that we're processing a FMD detector
340	Int_t detId = reader->GetDetectorID();
341	if (detId != (AliFMD::kBaseDDL >> 8)) {
342	Error("ReadAdcs", "Detector ID %d != %d",
343	detId, (AliFMD::kBaseDDL >> 8));
344	break;
46501dfb	345	}
4347b38f	346	// Figure out what detector we're deling with
	347	oldDetector = detector;
	348	switch (ddl) {
	349	case 0: detector = 1; break;
	350	case 1: detector = 2; break;
	351	case 2: detector = 3; break;
	352	default:
	353	Error("ReadAdcs", "Unknown DDL 0x%x for FMD", ddl);
	354	return kFALSE;
46501dfb	355	}
4347b38f	356	AliDebug(10, Form("Reading ADCs for 0x%x - That is FMD%d",
	357	reader->GetEquipmentId(), detector));
	358	}
	359	counts.Reset(-1);
	360	offset = 0;
	361	}
	362
	363	counts[offset] = input.Count();
	364	offset++;
	365
	366	AliDebug(10, Form("ADC of FMD%d%c[%2d,%3d] += %d",
	367	detector, input.Ring(), input.Sector(),
	368	input.Strip(), input.Count()));
	369	oldDetector = detector;
	370	}
	371	return kTRUE;
	372	}
88cb7938	373
4347b38f	374	//____________________________________________________________________
	375	void
	376	AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
	377	{
	378	// Process a digit. Derived classes can overload this member
	379	// function to do stuff to the digit. However, it should write the
	380	// ADC count to the internal cache
	381	//
	382	// fAdcs(detector - 1, ring, sector, strip) = counts;
	383	//
	384	// In this implementation, we count the number of strips below
	385	// threshold. This we do to later choose what kind of
	386	// reconstruction algorithm we'd like to use.
	387	//
	388	UShort_t detector = digit->Detector();
	389	Char_t ring = digit->Ring();
	390	UShort_t sector = digit->Sector();
	391	UShort_t strip = digit->Strip();
	392
	393	UShort_t counts = SubtractPedestal(digit);
	394
	395	fAdcs(detector - 1, ring, sector, strip) = counts;
	396	if (counts < fThreshold) fEmptyStrips++;
	397	fTotalStrips++;
	398	}
1a42d4f2	399
4347b38f	400	//____________________________________________________________________
	401	UShort_t
	402	AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
	403	{
	404	// Member function to subtract the pedestal from a digit
	405	// This implementation does nothing, but a derived class could over
	406	// load this to subtract a pedestal that was given in a database or
	407	// something like that.
3d44ce66	408
4347b38f	409	Int_t counts =
	410	TMath::Max(Int_t(digit->Count1() - fPedestal - 3 * fPedestalWidth), 0);
	411	if (digit->Count2() >= 0)
	412	counts +=
	413	TMath::Max(Int_t(digit->Count2() - fPedestal - 3 * fPedestalWidth), 0);
	414	if (digit->Count3() >= 0)
	415	counts +=
	416	TMath::Max(Int_t(digit->Count3() - fPedestal - 3 * fPedestalWidth), 0);
	417
	418	if (counts < 0) counts = 0;
	419	return UShort_t(counts);
88cb7938	420	}
dc8af42e	421
4347b38f	422	//____________________________________________________________________
	423	void
	424	AliFMDReconstructor::ReconstructFromCache(Float_t zVertex) const
	425	{
	426	// Based on the information in the cache, do the reconstruction.
	427	Int_t nRecon = 0;
	428	// Loop over the detectors
	429	for (Int_t i = 1; i <= 3; i++) {
	430	AliFMDSubDetector* sub = 0;
	431	switch (i) {
	432	case 1: sub = fFMD->GetFMD1(); break;
	433	case 2: sub = fFMD->GetFMD2(); break;
	434	case 3: sub = fFMD->GetFMD3(); break;
	435	}
	436	if (!sub) continue;
	437
	438	// Loop over the rings in the detector
	439	for (Int_t j = 0; j < 2; j++) {
	440	Float_t rZ = 0;
	441	AliFMDRing* r = 0;
	442	switch (j) {
	443	case 0: r = sub->GetInner(); rZ = sub->GetInnerZ(); break;
	444	case 1: r = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
	445	}
	446	if (!r) continue;
	447
	448	// Calculate low/high theta and eta
	449	// FIXME: Is this right?
	450	Float_t realZ = zVertex + rZ;
	451	Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
	452	Float_t thetaIn = TMath::ATan2(r->GetLowR(), realZ);
	453	Float_t etaOut = - TMath::Log(TMath::Tan(thetaOut / 2));
	454	Float_t etaIn = - TMath::Log(TMath::Tan(thetaIn / 2));
	455	if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
	456	Float_t tmp = etaIn;
	457	etaIn = etaOut;
	458	etaOut = tmp;
	459	}
	460
	461	//-------------------------------------------------------------
	462	//
	463	// Here starts poisson method
	464	//
	465	// Calculate eta step per strip, number of eta steps, number of
	466	// phi steps, and check the sign of the eta increment
	467	Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
	468	Int_t nEta = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
	469	Int_t nPhi = Int_t(360. / fDeltaPhi);
	470	Float_t sign = TMath::Sign(Float_t(1.), etaIn);
121a60bd	471
4347b38f	472	AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
	473	sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
	474
	475	// Loop over relevant phi values
	476	for (Int_t p = 0; p < nPhi; p++) {
	477	Float_t minPhi = p * fDeltaPhi;
	478	Float_t maxPhi = minPhi + fDeltaPhi;
	479	UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
	480	UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
	481
	482	AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
	483	minPhi, maxPhi, minSector, maxSector));
	484	// Loop over relevant eta values
	485	for (Int_t e = nEta; e >= 0; --e) {
	486	Float_t maxEta = etaIn - sign * e * fDeltaEta;
	487	Float_t minEta = maxEta - sign * fDeltaEta;
	488	if (sign > 0) minEta = TMath::Max(minEta, etaOut);
	489	else minEta = TMath::Min(minEta, etaOut);
	490	Float_t theta1 = 2 * TMath::ATan(TMath::Exp(-minEta));
	491	Float_t theta2 = 2 * TMath::ATan(TMath::Exp(-maxEta));
	492	Float_t minR = TMath::Abs(realZ * TMath::Tan(theta2));
	493	Float_t maxR = TMath::Abs(realZ * TMath::Tan(theta1));
	494	UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
	495	UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
	496
	497	AliDebug(10, Form(" Now in eta range %f, %f (strips %d, %d)\n"
	498	" [radii %f, %f, thetas %f, %f, sign %d]",
	499	minEta, maxEta, minStrip, maxStrip,
	500	minR, maxR, theta1, theta2, sign));
	501
	502	// Count number of empty strips
	503	Int_t emptyStrips = 0;
	504	for (Int_t sector = minSector; sector < maxSector; sector++)
	505	for (Int_t strip = minStrip; strip < maxStrip; strip++)
	506	if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip)
	507	< fThreshold) emptyStrips++;
	508
	509	// The total number of strips
	510	Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
	511
	512	// Log ratio of empty to total number of strips
	513	AliDebug(10, Form("Lambda= %d / %d = %f",
	514	emptyStrips, nTotal,
	515	Float_t(emptyStrips) / nTotal));
	516
	517	Double_t lambda = (emptyStrips > 0 ?
	518	- TMath::Log(Double_t(emptyStrips) / nTotal) :
	519	1);
	520
	521	// The reconstructed number of particles is then given by
	522	Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
	523
	524	// Add a AliFMDParticles to the reconstruction tree.
	525	new((*fParticles)[nRecon])
	526	AliFMDParticles(sub->GetId(), r->GetId(),
	527	minSector, maxSector, minStrip, maxStrip,
	528	minEta, maxEta, minPhi, maxPhi,
	529	reconstructed, AliFMDParticles::kPoission);
	530	nRecon++;
	531	} // phi
	532	} // eta
	533	} // ring
	534	} // detector
	535	}
536
537
538	//____________________________________________________________________
539	void
540	AliFMDReconstructor::FillESD(AliRunLoader* /fRunLoader/,
541	AliESD* /esd/) const
121a60bd	542	{
	543	// nothing to be done
	544
	545	}
	546