[u/mrichter/AliRoot.git] / FMD / AliFMDReconstruction.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 TTask that constructs ReconstParticles (reconstructed particles) 
// out of Digits
// 
//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
//////////////////////////////////////////////////////////////////////////////

// --- ROOT system ---
#include "TTask.h"
#include "TTree.h"
#include "TSystem.h"
#include "TFile.h"
#include "TROOT.h"
#include "TFolder.h"

// --- Standard library ---
#include <stdlib.h>
#include "Riostream.h"

// --- AliRoot header files ---

#include "AliFMDdigit.h"
#include "AliFMDhit.h"
#include "AliFMDReconstParticles.h"
#include "AliFMD.h"
#include "AliFMDv1.h"
#include "AliFMDReconstruction.h"
#include "AliRun.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
ClassImp(AliFMDReconstruction)

        
//____________________________________________________________________________ 

AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","") 
{
  fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
  // add Task to //root/Tasks folder
  TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
  roottasks->Add(this) ; 
}
//____________________________________________________________________________ 

AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstParticlesFile):TTask("AliFMDReconstruction","")
{
  fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
  fReconstParticlesFile=ReconstParticlesFile ;
  fHeadersFile=HeaderFile ;
  //add Task to //root/Tasks folder
  TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
  roottasks->Add(this) ;     
}

//____________________________________________________________________________ 

AliFMDReconstruction::~AliFMDReconstruction()
{
}
//----------------------------------------------------------------------------

void AliFMDReconstruction::Exec(Option_t *option) 
{ 
  /*
    Reconstruct nember of particles 
    in given group of pads for given FMDvolume
    determine by numberOfVolume , 
    numberOfMinSector,numberOfMaxSector,
    numberOfMinRing, numberOgMaxRing
    Reconstruction method choose dependence on number of empty pads  
  */


  cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
  Int_t const knumVolumes=5;
  Int_t const knumSectors=40;
  Int_t const knumRings=768;
  Int_t padADC[10][50][800];
  Float_t eta, etain,etaout,rad,theta;
  Int_t ivol, iSector, iRing;
  Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
  Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
  Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
  Int_t numberOfRings[5]={768,384,768,384,768};
  Int_t numberOfSectors[5]=  {20,40,20,40,20};
  Int_t numberOfEtaIntervals[5];
  // number of ring for boundary 0.1 eta
  TBranch *brDigits=0;

  AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
  TClonesArray *fReconParticles=fFMD->ReconParticles();
  if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
  cout<<" fNevents "<<fNevents<<endl;
  for(Int_t ievent=0;ievent<fNevents;ievent++)
    { 
      cout<<" ievent "<<ievent<<endl;
      for (Int_t i=0; i<knumVolumes; i++)
	for(Int_t j=0; j<knumSectors; j++)
	  for(Int_t ij=0; ij<knumRings; ij++)
	    padADC[i][j][ij]=0;                    //zhachem ???
      gAlice->GetEvent(ievent) ;
      if(gAlice->TreeH()==0) return; 
      if(gAlice->TreeD()==0) return;
	brDigits=gAlice->TreeD()->GetBranch("FMD");
	if (!brDigits){
	  cerr<<"EXEC Branch FMD digits not found"<<endl;
        return;
      } 
 
      if(gAlice->TreeR()==0) gAlice->MakeTree("R");
      //Make branches
      fFMD->MakeBranch("R");

      
      Int_t zeroADC=1;
 
      AliFMDdigit  *fmdDigit;
       if (fFMD)
	{
	  gAlice->TreeD()->GetEvent(0); 
	  TClonesArray *fFMDdigits=fFMD->Digits();
	  Int_t nDigits=fFMDdigits->GetEntries();
	  cout<<" nDigits "<<nDigits<<endl;
	   Int_t recParticles[6];
	   Int_t nRecPart=0 ;
	   Int_t zeroPads=0;
	   Int_t numberOfPads=0; //To avoid warning
	   Int_t pedestal;
	   Float_t channelWidth=(22400*50)/1024;
	   for (Int_t digit=0;digit<nDigits;digit++) 
	     {
	       fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);    
	       ivol=fmdDigit->Volume();
	       iSector=fmdDigit->NumberOfSector();
	       iRing=fmdDigit->NumberOfRing();
	       pedestal=Int_t(gRandom->Gaus(500,250));
	       padADC[ivol-1][iSector-1][iRing-1]=
		 fmdDigit->ADCsignal()
		 -Int_t(Float_t(pedestal)/channelWidth);
	       if (padADC[ivol-1][iSector-1][iRing-1]<0) 
		 padADC[ivol-1][iSector-1][iRing-1]=0;
	     } //digit loop
	   Int_t rmin=0; Int_t rmax=0; //To avoid warning
	   Int_t smin=0; Int_t smax=0; //To avoid warning
	   AliHeader *header = gAlice->GetHeader();
	   AliGenEventHeader* genHeader = header->GenEventHeader();
	   TArrayF *o = new TArrayF(3); 
	   genHeader->PrimaryVertex(*o);
	   Float_t zVertex=o->At(2);
 	   for (ivol=0; ivol<knumVolumes; ivol++)
	     {
	       Float_t realZ=zVertex+z[ivol];
	       theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
	       etain = - TMath::Log( TMath::Tan(theta/2.));
	       theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
	       etaout=- TMath::Log( TMath::Tan(theta/2.));
	       numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
	       eta=etain;
	       for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++) 
		 {
		   theta = 2.*TMath::ATan (TMath::Exp (-eta));
		   Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
		   rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
		   eta=eta-0.1;
		   theta = 2.*TMath::ATan (TMath::Exp (-eta));
		   rad = TMath::Abs(realZ) * (TMath::Tan (theta));
		   rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
		   zeroPads=0;
		   smin=0;
		   smax=numberOfSectors[ivol]; 
		   for (Int_t iring=rmin; iring<rmax; iring++) 
		     {
		       numberOfPads=(rmax-rmin)*(smax-smin);
		       for 
			 (Int_t isector=1;
			  isector<=numberOfSectors[ivol]; 
			  isector++) 
			 if(padADC[ivol][isector-1][iring-1]<zeroADC)
			   zeroPads++;
		     } //ring
		   Float_t zerosRatio= 
		     (Float_t)zeroPads/(Float_t)numberOfPads;
		   recParticles[0]=ivol;
		   recParticles[1]=smin;
		   recParticles[2]=smax;
		   recParticles[3]=rmin;
		   recParticles[4]=rmax;
		   if (zerosRatio>0.1 )
		     recParticles[5]=
		       DeterminationByPoisson
		       (padADC,ivol+1,rmin,rmax,smin,smax);
		   else
		     recParticles[5]=
		       DeterminationByThresholds
		       (padADC,ivol+1,rmin,rmax,smin,smax);
		   new((*fReconParticles)[nRecPart++]) 
		 AliFMDReconstParticles(recParticles);   	       
		 } // eta
	     } // volume
	   
	}//if FMD
       gAlice->TreeR()->Reset();
       gAlice->TreeR()->Fill(); 
       gAlice->TreeR()->Write(0,TObject::kOverwrite);
    } //event loop
  cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
};

//------------------------------------------------------------
Int_t AliFMDReconstruction::DeterminationByThresholds
(Int_t padAdc[][50][800],Int_t volume, Int_t rmin, Int_t rmax, 
 Int_t smin, Int_t smax)
{
  /*
    reconstruct number of particles by 
    energy deposited threshold method 
    Using if number of empty pads less then 10%

*/
  cout<<"\nStart threshold method\n";
  
  Int_t thresholdInner[30]={
    0,     14,  28,    42,   57,     
    72,    89,  104,  124,  129, 
    164,  174,  179,  208,  228, 
    248,  268,   317,  337,  357, 
    392,  407,  416,  426,  436, 
    461,  468,  493,  506,  515}; 

  Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
			    132, 165, 198, 231, 
			    264, 286, 308, 334, 
			    352, 374, 418, 440,
			    462, 484, 506, 528,
			    550, 572, 594, 616}; 
  Int_t threshold[30];
  for (Int_t it=0; it<30; it++) {
    if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
    if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
  }
  Int_t thresholdArraySize = 30;     
  Int_t numPart=0;
  //Inner Si
  for (Int_t iring=rmin; iring<rmax; iring++) 
    {
      for (Int_t isector=smin; isector<smax; isector++) 
	{
	  for (int i=0;i<thresholdArraySize-1;i++)
	    {	      
	      if(padAdc[volume-1][isector][iring]>threshold[i]
		 &&padAdc[volume-1][isector][iring]<=threshold[i+1])
		{
		  numPart+=i;
		  break;
		}; // if in threshol interval
	    }; //threshold_array_size
	}; //iring
    }; //sector
  cout<<"\nEnd threshol method"<<endl;
  return numPart;
}


 //__________________________________________________________________

Int_t AliFMDReconstruction::DeterminationByPoisson 
(Int_t padAdc[][50][800],
 Int_t vol, Int_t rmin, Int_t rmax, 
 Int_t secmin, Int_t secmax)
{
  /* 
     reconstruct number of particles by Poisson statistic method 
     according number of empty pad in chosen group of pads 
     using if number of empty pads more then 10% 

  */

  //  cout<<"\nStart Poisson method";
  Int_t thresholdADC=1;   
  Int_t zeropads=0;
  for (Int_t i=rmin;i<rmax;i++)
    {
      for (Int_t j=secmin;j<secmax;j++)
	{
	  if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
	};
    };
  Float_t lambda=-TMath::Log(Float_t(zeropads)/
			     ( Float_t(secmax-secmin)*
			       Float_t(rmax-rmin))); //+1 zdes' ne nado
  Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
  //  cout<<"\nEnd Poisson method"<<endl;
  return fRecon;
};

//------------------------------------------------------------------
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
	16	//_________________________________________________________________________
	17	// This is a TTask that constructs ReconstParticles (reconstructed particles)
	18	// out of Digits
	19	//
	20	//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
	21	//////////////////////////////////////////////////////////////////////////////
	22
dc8af42e	23	// --- ROOT system ---
	24	#include "TTask.h"
	25	#include "TTree.h"
	26	#include "TSystem.h"
	27	#include "TFile.h"
37c55dc0	28	#include "TROOT.h"
	29	#include "TFolder.h"
	30
dc8af42e	31	// --- Standard library ---
37c55dc0	32	#include <stdlib.h>
e75fab8f	33	#include "Riostream.h"
dc8af42e	34
	35	// --- AliRoot header files ---
	36
	37	#include "AliFMDdigit.h"
1a42d4f2	38	#include "AliFMDhit.h"
dc8af42e	39	#include "AliFMDReconstParticles.h"
	40	#include "AliFMD.h"
	41	#include "AliFMDv1.h"
	42	#include "AliFMDReconstruction.h"
	43	#include "AliRun.h"
1a42d4f2	44	#include "AliHeader.h"
1a42d4f2	45	#include "AliGenEventHeader.h"
dc8af42e	46	ClassImp(AliFMDReconstruction)
	47
	48
	49	//____________________________________________________________________________
	50
	51	AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","")
	52	{
	53	fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
	54	// add Task to //root/Tasks folder
	55	TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
	56	roottasks->Add(this) ;
	57	}
	58	//____________________________________________________________________________
	59
	60	AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstParticlesFile):TTask("AliFMDReconstruction","")
	61	{
	62	fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
	63	fReconstParticlesFile=ReconstParticlesFile ;
	64	fHeadersFile=HeaderFile ;
	65	//add Task to //root/Tasks folder
	66	TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
	67	roottasks->Add(this) ;
	68	}
	69
	70	//____________________________________________________________________________
	71
	72	AliFMDReconstruction::~AliFMDReconstruction()
	73	{
dc8af42e	74	}
cb1df35e	75	//----------------------------------------------------------------------------
dc8af42e	76
37c55dc0	77	void AliFMDReconstruction::Exec(Option_t *option)
dc8af42e	78	{
a39f3926	79	/*
	80	Reconstruct nember of particles
	81	in given group of pads for given FMDvolume
	82	determine by numberOfVolume ,
	83	numberOfMinSector,numberOfMaxSector,
	84	numberOfMinRing, numberOgMaxRing
	85	Reconstruction method choose dependence on number of empty pads
	86	*/
	87
	88
1a42d4f2	89	cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
a39f3926	90	Int_t const knumVolumes=5;
	91	Int_t const knumSectors=40;
	92	Int_t const knumRings=768;
	93	Int_t padADC[10][50][800];
1a42d4f2	94	Float_t eta, etain,etaout,rad,theta;
cb1df35e	95	Int_t ivol, iSector, iRing;
1a42d4f2	96	Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
	97	Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
	98	Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
a39f3926	99	Int_t numberOfRings[5]={768,384,768,384,768};
	100	Int_t numberOfSectors[5]= {20,40,20,40,20};
	101	Int_t numberOfEtaIntervals[5];
cb1df35e	102	// number of ring for boundary 0.1 eta
1a42d4f2	103	TBranch *brDigits=0;
cb1df35e	104
a39f3926	105	AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
a39f3926	106	TClonesArray *fReconParticles=fFMD->ReconParticles();
cb1df35e	107	if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
1a42d4f2	108	cout<<" fNevents "<<fNevents<<endl;
dc8af42e	109	for(Int_t ievent=0;ievent<fNevents;ievent++)
dc8af42e	110	{
1a42d4f2	111	cout<<" ievent "<<ievent<<endl;
a39f3926	112	for (Int_t i=0; i<knumVolumes; i++)
	113	for(Int_t j=0; j<knumSectors; j++)
	114	for(Int_t ij=0; ij<knumRings; ij++)
	115	padADC[i][j][ij]=0; //zhachem ???
dc8af42e	116	gAlice->GetEvent(ievent) ;
dc8af42e	117	if(gAlice->TreeH()==0) return;
1a42d4f2	118	if(gAlice->TreeD()==0) return;
	119	brDigits=gAlice->TreeD()->GetBranch("FMD");
	120	if (!brDigits){
	121	cerr<<"EXEC Branch FMD digits not found"<<endl;
	122	return;
	123	}
	124
dc8af42e	125	if(gAlice->TreeR()==0) gAlice->MakeTree("R");
dc8af42e	126	//Make branches
a39f3926	127	fFMD->MakeBranch("R");
1a42d4f2	128
37c55dc0	129
cb1df35e	130	Int_t zeroADC=1;
1a42d4f2	131
1a42d4f2	132	AliFMDdigit *fmdDigit;
a39f3926	133	if (fFMD)
dc8af42e	134	{
dc8af42e	135	gAlice->TreeD()->GetEvent(0);
a39f3926	136	TClonesArray *fFMDdigits=fFMD->Digits();
a39f3926	137	Int_t nDigits=fFMDdigits->GetEntries();
1a42d4f2	138	cout<<" nDigits "<<nDigits<<endl;
a39f3926	139	Int_t recParticles[6];
dc8af42e	140	Int_t nRecPart=0 ;
a39f3926	141	Int_t zeroPads=0;
a39f3926	142	Int_t numberOfPads=0; //To avoid warning
cb1df35e	143	Int_t pedestal;
cb1df35e	144	Float_t channelWidth=(22400*50)/1024;
dc8af42e	145	for (Int_t digit=0;digit<nDigits;digit++)
dc8af42e	146	{
a39f3926	147	fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);
cb1df35e	148	ivol=fmdDigit->Volume();
	149	iSector=fmdDigit->NumberOfSector();
	150	iRing=fmdDigit->NumberOfRing();
	151	pedestal=Int_t(gRandom->Gaus(500,250));
a39f3926	152	padADC[ivol-1][iSector-1][iRing-1]=
cb1df35e	153	fmdDigit->ADCsignal()
cb1df35e	154	-Int_t(Float_t(pedestal)/channelWidth);
a39f3926	155	if (padADC[ivol-1][iSector-1][iRing-1]<0)
a39f3926	156	padADC[ivol-1][iSector-1][iRing-1]=0;
dc8af42e	157	} //digit loop
a39f3926	158	Int_t rmin=0; Int_t rmax=0; //To avoid warning
a39f3926	159	Int_t smin=0; Int_t smax=0; //To avoid warning
1a42d4f2	160	AliHeader *header = gAlice->GetHeader();
	161	AliGenEventHeader* genHeader = header->GenEventHeader();
	162	TArrayF *o = new TArrayF(3);
	163	genHeader->PrimaryVertex(*o);
	164	Float_t zVertex=o->At(2);
a39f3926	165	for (ivol=0; ivol<knumVolumes; ivol++)
cb1df35e	166	{
1a42d4f2	167	Float_t realZ=zVertex+z[ivol];
	168	theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
	169	etain = - TMath::Log( TMath::Tan(theta/2.));
	170	theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
	171	etaout=- TMath::Log( TMath::Tan(theta/2.));
a39f3926	172	numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
1a42d4f2	173	eta=etain;
a39f3926	174	for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
cb1df35e	175	{
1a42d4f2	176	theta = 2.*TMath::ATan (TMath::Exp (-eta));
1a42d4f2	177	Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
a39f3926	178	rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
1a42d4f2	179	eta=eta-0.1;
	180	theta = 2.*TMath::ATan (TMath::Exp (-eta));
	181	rad = TMath::Abs(realZ) * (TMath::Tan (theta));
a39f3926	182	rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
	183	zeroPads=0;
	184	smin=0;
	185	smax=numberOfSectors[ivol];
	186	for (Int_t iring=rmin; iring<rmax; iring++)
cb1df35e	187	{
a39f3926	188	numberOfPads=(rmax-rmin)*(smax-smin);
cb1df35e	189	for
cb1df35e	190	(Int_t isector=1;
a39f3926	191	isector<=numberOfSectors[ivol];
cb1df35e	192	isector++)
a39f3926	193	if(padADC[ivol][isector-1][iring-1]<zeroADC)
a39f3926	194	zeroPads++;
cb1df35e	195	} //ring
cb1df35e	196	Float_t zerosRatio=
a39f3926	197	(Float_t)zeroPads/(Float_t)numberOfPads;
	198	recParticles[0]=ivol;
	199	recParticles[1]=smin;
	200	recParticles[2]=smax;
	201	recParticles[3]=rmin;
	202	recParticles[4]=rmax;
	203	if (zerosRatio>0.1 )
	204	recParticles[5]=
	205	DeterminationByPoisson
	206	(padADC,ivol+1,rmin,rmax,smin,smax);
cb1df35e	207	else
a39f3926	208	recParticles[5]=
	209	DeterminationByThresholds
	210	(padADC,ivol+1,rmin,rmax,smin,smax);
cb1df35e	211	new((*fReconParticles)[nRecPart++])
a39f3926	212	AliFMDReconstParticles(recParticles);
1a42d4f2	213	} // eta
1a42d4f2	214	} // volume
cb1df35e	215
cb1df35e	216	}//if FMD
1a42d4f2	217	gAlice->TreeR()->Reset();
	218	gAlice->TreeR()->Fill();
	219	gAlice->TreeR()->Write(0,TObject::kOverwrite);
dc8af42e	220	} //event loop
1a42d4f2	221	cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
cb1df35e	222	};
dc8af42e	223
cb1df35e	224	//------------------------------------------------------------
a39f3926	225	Int_t AliFMDReconstruction::DeterminationByThresholds
	226	(Int_t padAdc[][50][800],Int_t volume, Int_t rmin, Int_t rmax,
	227	Int_t smin, Int_t smax)
dc8af42e	228	{
a39f3926	229	/*
	230	reconstruct number of particles by
	231	energy deposited threshold method
	232	Using if number of empty pads less then 10%
	233
	234	*/
1a42d4f2	235	cout<<"\nStart threshold method\n";
cb1df35e	236
	237	Int_t thresholdInner[30]={
	238	0, 14, 28, 42, 57,
	239	72, 89, 104, 124, 129,
	240	164, 174, 179, 208, 228,
	241	248, 268, 317, 337, 357,
	242	392, 407, 416, 426, 436,
	243	461, 468, 493, 506, 515};
	244
	245	Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
	246	132, 165, 198, 231,
	247	264, 286, 308, 334,
	248	352, 374, 418, 440,
	249	462, 484, 506, 528,
	250	550, 572, 594, 616};
	251	Int_t threshold[30];
	252	for (Int_t it=0; it<30; it++) {
	253	if(volume==1\|\|volume==3\|\|volume==5) threshold[it]=thresholdInner[it];
	254	if(volume==2\|\|volume==4) threshold[it]=thresholdOuter[it];
	255	}
a39f3926	256	Int_t thresholdArraySize = 30;
a39f3926	257	Int_t numPart=0;
cb1df35e	258	//Inner Si
a39f3926	259	for (Int_t iring=rmin; iring<rmax; iring++)
cb1df35e	260	{
a39f3926	261	for (Int_t isector=smin; isector<smax; isector++)
cb1df35e	262	{
a39f3926	263	for (int i=0;i<thresholdArraySize-1;i++)
cb1df35e	264	{
a39f3926	265	if(padAdc[volume-1][isector][iring]>threshold[i]
a39f3926	266	&&padAdc[volume-1][isector][iring]<=threshold[i+1])
cb1df35e	267	{
a39f3926	268	numPart+=i;
cb1df35e	269	break;
	270	}; // if in threshol interval
	271	}; //threshold_array_size
	272	}; //iring
	273	}; //sector
1a42d4f2	274	cout<<"\nEnd threshol method"<<endl;
a39f3926	275	return numPart;
dc8af42e	276	}
cb1df35e	277
	278
	279	//__________________________________________________________________
	280
a39f3926	281	Int_t AliFMDReconstruction::DeterminationByPoisson
	282	(Int_t padAdc[][50][800],
	283	Int_t vol, Int_t rmin, Int_t rmax,
	284	Int_t secmin, Int_t secmax)
dc8af42e	285	{
a39f3926	286	/*
	287	reconstruct number of particles by Poisson statistic method
	288	according number of empty pad in chosen group of pads
	289	using if number of empty pads more then 10%
	290
	291	*/
	292
1a42d4f2	293	// cout<<"\nStart Poisson method";
a39f3926	294	Int_t thresholdADC=1;
cb1df35e	295	Int_t zeropads=0;
	296	for (Int_t i=rmin;i<rmax;i++)
	297	{
	298	for (Int_t j=secmin;j<secmax;j++)
	299	{
a39f3926	300	if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
cb1df35e	301	};
	302	};
	303	Float_t lambda=-TMath::Log(Float_t(zeropads)/
	304	( Float_t(secmax-secmin)*
	305	Float_t(rmax-rmin))); //+1 zdes' ne nado
a39f3926	306	Int_t fRecon=(Int_t)(lambda(secmax-secmin)(rmax-rmin)); //+1 zdes' ne nado
cb1df35e	307	// cout<<"\nEnd Poisson method"<<endl;
a39f3926	308	return fRecon;
cb1df35e	309	};
	310
	311	//------------------------------------------------------------------
	312
	313
	314
	315
	316
dc8af42e	317
	318
	319
37c55dc0	320
	321
	322
	323
	324
	325
	326
	327
	328
	329
	330
	331
	332
	333
	334
	335
	336
	337
1a42d4f2	338
dc8af42e	339
dc8af42e	340