[u/mrichter/AliRoot.git] / FMD / AliFMDMultPoisson.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$ */

//____________________________________________________________________
// 
// Reconstruct charged particle multiplicity in the FMD 
// 
// [See also the AliFMDReconstructor class]
// 
// This class reconstructs the muliplicity in regions based on the
// ratio of empty to full strips. 
// Hence the name `poisson'
//
#include "AliFMD.h"			// ALIFMD_H
#include "AliFMDGeometry.h"		// ALIFMDGEOMETRY_H
#include "AliFMDDetector.h"		// ALIFMDDETECTOR_H
#include "AliFMDRing.h"			// ALIFMDRING_H
#include "AliFMDMultPoisson.h"		// ALIFMDMULTPOISSON_H
#include "AliFMDMultRegion.h"		// ALIFMDMULTREGION_H
#include "AliFMDDigit.h"		// ALIFMDDIGIT_H
#include "AliLog.h"			// ALILOG_H
#include <TClonesArray.h>               // ROOT_TClonesArray
#include <TTree.h>               	// ROOT_TTree

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

//____________________________________________________________________
AliFMDMultPoisson::AliFMDMultPoisson()
  : AliFMDMultAlgorithm("Poisson", "Poisson"),
    fDeltaEta(0), 
    fDeltaPhi(0), 
    fThreshold(0)
{
  // CTOR
  SetDeltaEta();
  SetDeltaPhi();
  SetThreshold();
  fMult = new TClonesArray("AliFMDMultRegion", 1000);
}

//____________________________________________________________________
void
AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ) 
{
  // Reset internal data
  AliFMDMultAlgorithm::PreEvent(tree, ipZ);
  fCurrentVertexZ = ipZ;
  fEmpty.Reset(kFALSE);

  // Make a branch in the reconstruction tree. 
  const Int_t kBufferSize = 16000;
  fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);  
  
}

//____________________________________________________________________
void
AliFMDMultPoisson::ProcessDigit(AliFMDDigit*  digit, 
				Float_t       /* eta */, 
				Float_t       /* phi */, 
				UShort_t      count)
{
  // Process one digit. 
  // 
  // Parameters: 
  //    
  //   digit		Digit to process 
  //   ipZ		Z--coordinate of the primary interaction
  //                    vertex of this event 
  //
  if (!digit) return;
  if (count < fThreshold) fEmpty(digit->Detector() - 1, 
				 digit->Ring(), 
				 digit->Sector(), 
				 digit->Strip()) = kTRUE;
}

//____________________________________________________________________
void
AliFMDMultPoisson::PostEvent() 
{
  // Fill the branch 
  // Based on the information in the cache, do the reconstruction. 

  // Loop over the detectors 
  for (Int_t i = 1; i <= 3; i++) {
    AliFMDGeometry* fmd = AliFMDGeometry::Instance();
    AliFMDDetector* sub = fmd->GetDetector(i);
    if (!sub) continue;
	
    // Loop over the rings in the detector
    for (Int_t j = 0; j < 2; j++) {
      AliFMDRing* r  = sub->GetRing((j == 0 ? 'I' : 'O'));
      Float_t     rZ = sub->GetRingZ((j == 0 ? 'I' : 'O'));
      if (!r) continue;
      
      // Calculate low/high theta and eta 
      // FIXME: Is this right? 
      Float_t realZ    = fCurrentVertexZ + 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));
	  // Calculate the weighted mean eta of the region
	  Float_t  minW2    = TMath::Power(minR * 2 * TMath::Pi() * 
					   ((maxPhi - minPhi)/360),2);
	  Float_t  maxW2    = TMath::Power(minR * 2 * TMath::Pi() * 
					   ((maxPhi - minPhi)/360), 2);
	  Float_t  meanEta  = ((minEta / minW2 + maxEta / maxW2) / 
			       (1 / (minW2 + maxW2)));
	  //UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
	  // UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);

	  UShort_t minStrip = UShort_t(r->GetNStrips() - 
				       (etaIn - minEta) * stripEta + 0.5);
	  UShort_t maxStrip = UShort_t(r->GetNStrips() - 
				       (etaIn - maxEta) * 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 (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip)) 
		emptyStrips++;
	  
	  // The total number of strips 
	  Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
	  // Log ratio of empty to total number of strips 
	  
	  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);
	  AliDebug(10, Form("Lambda= %d / %f = %f particles %d", 
			    emptyStrips, nTotal, 
			    Float_t(emptyStrips) / nTotal , reconstructed));
	    
	  // Add a AliFMDMultRegion to the reconstruction tree. 
	  AliFMDMultRegion* m = new((*fMult)[fNMult])   
	    AliFMDMultRegion(sub->GetId(), r->GetId(),
			     minSector, maxSector, minStrip, maxStrip,
			     minEta, maxEta, meanEta, minPhi, maxPhi,
			     reconstructed, AliFMDMultRegion::kPoission);
	  (void)m;
	  fNMult++;
	} // phi 
      } // eta
    } // ring 
  } // detector 
}


//____________________________________________________________________
// 
// EOF
//
Commit	Line	Data
56b1929b	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
	16	/* $Id$ */
	17
	18	//____________________________________________________________________
	19	//
7684b53c	20	// Reconstruct charged particle multiplicity in the FMD
	21	//
	22	// [See also the AliFMDReconstructor class]
56b1929b	23	//
7684b53c	24	// This class reconstructs the muliplicity in regions based on the
7684b53c	25	// ratio of empty to full strips.
088f8e79	26	// Hence the name `poisson'
7684b53c	27	//
56b1929b	28	#include "AliFMD.h" // ALIFMD_H
1a1fdef7	29	#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
	30	#include "AliFMDDetector.h" // ALIFMDDETECTOR_H
	31	#include "AliFMDRing.h" // ALIFMDRING_H
56b1929b	32	#include "AliFMDMultPoisson.h" // ALIFMDMULTPOISSON_H
	33	#include "AliFMDMultRegion.h" // ALIFMDMULTREGION_H
	34	#include "AliFMDDigit.h" // ALIFMDDIGIT_H
	35	#include "AliLog.h" // ALILOG_H
	36	#include <TClonesArray.h> // ROOT_TClonesArray
	37	#include <TTree.h> // ROOT_TTree
	38
	39	//____________________________________________________________________
925e6570	40	ClassImp(AliFMDMultPoisson)
1a1fdef7	41	#if 0
	42	; // This is here to keep Emacs for indenting the next line
	43	#endif
56b1929b	44
	45	//____________________________________________________________________
	46	AliFMDMultPoisson::AliFMDMultPoisson()
	47	: AliFMDMultAlgorithm("Poisson", "Poisson"),
	48	fDeltaEta(0),
	49	fDeltaPhi(0),
	50	fThreshold(0)
	51	{
088f8e79	52	// CTOR
56b1929b	53	SetDeltaEta();
	54	SetDeltaPhi();
	55	SetThreshold();
	56	fMult = new TClonesArray("AliFMDMultRegion", 1000);
	57	}
	58
	59	//____________________________________________________________________
	60	void
	61	AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ)
	62	{
	63	// Reset internal data
	64	AliFMDMultAlgorithm::PreEvent(tree, ipZ);
	65	fCurrentVertexZ = ipZ;
69b696b9	66	fEmpty.Reset(kFALSE);
56b1929b	67
	68	// Make a branch in the reconstruction tree.
	69	const Int_t kBufferSize = 16000;
	70	fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);
	71
	72	}
	73
	74	//____________________________________________________________________
	75	void
	76	AliFMDMultPoisson::ProcessDigit(AliFMDDigit* digit,
	77	Float_t /* eta */,
	78	Float_t /* phi */,
	79	UShort_t count)
	80	{
	81	// Process one digit.
	82	//
	83	// Parameters:
	84	//
	85	// digit Digit to process
	86	// ipZ Z--coordinate of the primary interaction
	87	// vertex of this event
	88	//
	89	if (!digit) return;
	90	if (count < fThreshold) fEmpty(digit->Detector() - 1,
	91	digit->Ring(),
	92	digit->Sector(),
	93	digit->Strip()) = kTRUE;
	94	}
	95
	96	//____________________________________________________________________
	97	void
	98	AliFMDMultPoisson::PostEvent()
	99	{
	100	// Fill the branch
	101	// Based on the information in the cache, do the reconstruction.
	102
	103	// Loop over the detectors
	104	for (Int_t i = 1; i <= 3; i++) {
1a1fdef7	105	AliFMDGeometry* fmd = AliFMDGeometry::Instance();
1a1fdef7	106	AliFMDDetector* sub = fmd->GetDetector(i);
56b1929b	107	if (!sub) continue;
	108
	109	// Loop over the rings in the detector
	110	for (Int_t j = 0; j < 2; j++) {
1a1fdef7	111	AliFMDRing* r = sub->GetRing((j == 0 ? 'I' : 'O'));
1a1fdef7	112	Float_t rZ = sub->GetRingZ((j == 0 ? 'I' : 'O'));
56b1929b	113	if (!r) continue;
	114
	115	// Calculate low/high theta and eta
	116	// FIXME: Is this right?
	117	Float_t realZ = fCurrentVertexZ + rZ;
	118	Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
	119	Float_t thetaIn = TMath::ATan2(r->GetLowR(), realZ);
	120	Float_t etaOut = - TMath::Log(TMath::Tan(thetaOut / 2));
	121	Float_t etaIn = - TMath::Log(TMath::Tan(thetaIn / 2));
	122	if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
	123	Float_t tmp = etaIn;
	124	etaIn = etaOut;
	125	etaOut = tmp;
	126	}
	127
	128	//-------------------------------------------------------------
	129	//
	130	// Here starts poisson method
	131	//
	132	// Calculate eta step per strip, number of eta steps, number of
	133	// phi steps, and check the sign of the eta increment
	134	Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
	135	Int_t nEta = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
	136	Int_t nPhi = Int_t(360. / fDeltaPhi);
	137	Float_t sign = TMath::Sign(Float_t(1.), etaIn);
	138
	139	AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
	140	sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
	141
	142	// Loop over relevant phi values
	143	for (Int_t p = 0; p < nPhi; p++) {
	144	Float_t minPhi = p * fDeltaPhi;
	145	Float_t maxPhi = minPhi + fDeltaPhi;
	146	UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
	147	UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
	148
c263f0ab	149	// AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
c263f0ab	150	// minPhi, maxPhi, minSector, maxSector));
56b1929b	151	// Loop over relevant eta values
	152	for (Int_t e = nEta; e >= 0; --e) {
	153	Float_t maxEta = etaIn - sign * e * fDeltaEta;
	154	Float_t minEta = maxEta - sign * fDeltaEta;
	155	if (sign > 0) minEta = TMath::Max(minEta, etaOut);
	156	else minEta = TMath::Min(minEta, etaOut);
	157	Float_t theta1 = 2 * TMath::ATan(TMath::Exp(-minEta));
	158	Float_t theta2 = 2 * TMath::ATan(TMath::Exp(-maxEta));
	159	Float_t minR = TMath::Abs(realZ * TMath::Tan(theta2));
	160	Float_t maxR = TMath::Abs(realZ * TMath::Tan(theta1));
a3537838	161	// Calculate the weighted mean eta of the region
	162	Float_t minW2 = TMath::Power(minR * 2 * TMath::Pi() *
	163	((maxPhi - minPhi)/360),2);
	164	Float_t maxW2 = TMath::Power(minR * 2 * TMath::Pi() *
	165	((maxPhi - minPhi)/360), 2);
	166	Float_t meanEta = ((minEta / minW2 + maxEta / maxW2) /
	167	(1 / (minW2 + maxW2)));
c263f0ab	168	//UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
	169	// UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
	170
a3537838	171	UShort_t minStrip = UShort_t(r->GetNStrips() -
	172	(etaIn - minEta) * stripEta + 0.5);
	173	UShort_t maxStrip = UShort_t(r->GetNStrips() -
	174	(etaIn - maxEta) * stripEta + 0.5);
56b1929b	175
56b1929b	176	AliDebug(10, Form(" Now in eta range %f, %f (strips %d, %d)\n"
c263f0ab	177	" [radii %f, %f, thetas %f, %f, sign %d]",
c263f0ab	178	minEta, maxEta, minStrip, maxStrip,
56b1929b	179	minR, maxR, theta1, theta2, sign));
c263f0ab	180
56b1929b	181	// Count number of empty strips
	182	Int_t emptyStrips = 0;
	183	for (Int_t sector = minSector; sector < maxSector; sector++)
	184	for (Int_t strip = minStrip; strip < maxStrip; strip++)
	185	if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip))
	186	emptyStrips++;
	187
	188	// The total number of strips
	189	Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
56b1929b	190	// Log ratio of empty to total number of strips
56b1929b	191
	192	Double_t lambda = (emptyStrips > 0 ?
	193	- TMath::Log(Double_t(emptyStrips) / nTotal) :
	194	1);
	195
	196	// The reconstructed number of particles is then given by
	197	Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
c263f0ab	198	AliDebug(10, Form("Lambda= %d / %f = %f particles %d",
	199	emptyStrips, nTotal,
	200	Float_t(emptyStrips) / nTotal , reconstructed));
56b1929b	201
56b1929b	202	// Add a AliFMDMultRegion to the reconstruction tree.
7684b53c	203	AliFMDMultRegion* m = new((*fMult)[fNMult])
56b1929b	204	AliFMDMultRegion(sub->GetId(), r->GetId(),
56b1929b	205	minSector, maxSector, minStrip, maxStrip,
a3537838	206	minEta, maxEta, meanEta, minPhi, maxPhi,
56b1929b	207	reconstructed, AliFMDMultRegion::kPoission);
7684b53c	208	(void)m;
56b1929b	209	fNMult++;
	210	} // phi
	211	} // eta
	212	} // ring
	213	} // detector
	214	}
	215
	216
	217	//____________________________________________________________________
	218	//
	219	// EOF
	220	//