[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 "AliRunLoader.h"
#include "AliLoader.h"

#include "AliFMDdigit.h"
#include "AliFMDhit.h"
#include "AliFMDReconstParticles.h"
#include "AliFMD.h"
#include "AliFMDv1.h"
#include "AliFMDReconstruction.h"
#include "AliRun.h"
#include "AliConfig.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
  fRunLoader = 0x0;
  
}
//____________________________________________________________________________ 

AliFMDReconstruction::AliFMDReconstruction(AliRunLoader* rl):TTask("AliFMDReconstruction","")
{

  if (rl == 0x0)
   {
     Fatal("AliFMDReconstruction","Argument AliRunLoader* is null!");
     return;
   }
   
  fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file

  fRunLoader = rl;
  AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
  if (gime == 0x0)
   {
     Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
     return;//never reached
   }
  //add Task to //root/Tasks folder
  gime->PostReconstructioner(this);
}

//____________________________________________________________________________ 

AliFMDReconstruction::~AliFMDReconstruction()
{
}

//____________________________________________________________________________

void AliFMDReconstruction::Exec(Option_t *option) 
{ 
 //Collects all digits in the same active volume into number of particles
  /*
    Reconstruct number 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

  
  if (fRunLoader == 0x0)
   {
    Error("Exec","Run Loader loader is NULL - Session not opened");
    return;
   }
  AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
  if (gime == 0x0)
   {
     Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
     return;//never reached
   }
   
  fRunLoader->LoadgAlice();
  Int_t retval;
  
  retval = gime->LoadHits("READ"); 
  if (retval)
   {
     Error("Exec","Error occured while loading hits. Exiting.");
     return;
   }

  retval = gime->LoadDigits("READ"); 
  if (retval)
   {
     Error("Exec","Error occured while loading digits. Exiting.");
     return;
   }
  
  AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
  TClonesArray *fReconParticles=fFMD->ReconParticles();
  
  TTree* treeD = gime->TreeD();
  if (treeD == 0x0)
   {
     Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
     return;
   }
  if(fNevents == 0) fNevents=(Int_t)treeD->GetEntries();
  //PH Do we use TreeE (kinematics), or TreeD (digits) toaccess the number
  //PH of events?
//PH   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
//PH   cout<<" fNevents "<<fNevents<<endl;

  TBranch *brDigits=0;

  for(Int_t ievent=0;ievent<fNevents;ievent++)
    { 
      fRunLoader->GetEvent(ievent) ;

      TTree* treeH = gime->TreeH();
      if (treeH == 0x0)
       {
         Error("Exec","Can not get TreeH");
         return;
       }
      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 ???

	brDigits=treeD->GetBranch("FMD");
	if (!brDigits){
	  cerr<<"EXEC Branch FMD digits not found"<<endl;
        return;
      } 
 
      if(gime->TreeR()==0) gime->MakeTree("R");

      //Make branches
      fFMD->MakeBranch("R");

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