[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
// 
//-- 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"
#include "TH2F.h"
#include "TRandom.h"

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

// --- AliRoot header files ---

#include "AliRunLoader.h"
#include "AliLoader.h"

#include "AliFMDdigit.h"
#include "AliFMDReconstParticles.h"
#include "AliFMDReconstructor.h"
#include "AliRun.h"
#include "AliConfig.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"

ClassImp(AliFMDReconstructor)

        
//____________________________________________________________________________

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
    determine by numberOfVolume , 
    numberOfMinSector,numberOfMaxSector,
    numberOfMinRing, numberOgMaxRing
    Reconstructor method choose dependence on number of empty pads  
  */


#ifdef DEBUG
  Info("Exec ","Start");
#endif

  Int_t const knumVolumes=5;
  Int_t padADC;
  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]={512,256,512,256,512};
  Int_t numberOfSectors[5]=  {20,40,20,40,20};
  Int_t numberOfEtaIntervals[5];
  // number of ring for boundary 0.1 eta

  
  if (runLoader == 0x0)
   {
    Error("Exec","Run Loader loader is NULL - Session not opened");
    return;
   }

  AliLoader* plFMD = runLoader->GetLoader("FMDLoader");
  if (plFMD == 0x0)
   {
     Fatal("AliFMDReconstructor","Can not find FMD (loader) in specified event");
     return;//never reached
   }
   
  if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  if (!runLoader->TreeE()) runLoader->LoadHeader();

  TDirectory* cwd = gDirectory;
  gDirectory = 0x0;
  Text_t buf1[20];
  TH2F* hTotal[10];
  for (Int_t j=1; j<=5; j++){
    sprintf(buf1,"hTotal%d",j);
    
    hTotal[j] = new TH2F(buf1," Number of primary particles ",
			 numberOfSectors[j-1],1,numberOfSectors[j-1],
			 numberOfRings[j-1],1,numberOfRings[j-1]);
  }   
  gDirectory = cwd;
 
  plFMD->LoadRecPoints("RECREATE");
  TClonesArray* reconParticles = new TClonesArray("AliFMDReconstParticles"); 
  TClonesArray* digits = new TClonesArray("AliFMDdigit");
 
  Int_t retval=0;     
  Int_t nevents=Int_t (runLoader->TreeE()->GetEntries()); 
#ifdef DEBUG
  cout<<" nevents "<<nevents<<endl;
#endif
   for(Int_t ievent=0;ievent<nevents;ievent++)
    { 
#ifdef DEBUG
      cout<<" *** event "<<ievent<<endl; 
#endif
      runLoader->GetEvent(ievent) ;
      //event z-vertex for correction eta-rad dependence      
      AliHeader *header = runLoader->GetHeader();
      AliGenEventHeader* genHeader = header->GenEventHeader();
      TArrayF *o = new TArrayF(3); 
      if (genHeader) genHeader->PrimaryVertex(*o);
      Float_t zVertex=o->At(2);

      for (Int_t i=0; i<5; i++)
        hTotal[i+1]->Reset();
      
      retval = plFMD->LoadDigits("READ"); 
      if (retval)
	{
	  Error("Exec","Error occured while loading digits. Exiting.");
	  return;
	}
      
      TTree* treeD = plFMD->TreeD();
      if (treeD == 0x0)
	{
	  Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
	  return;
	}
      
      TBranch *brDigits=0;
	    
      brDigits=treeD->GetBranch("FMD");

      if (brDigits) {
	brDigits->SetAddress(&digits);
      }else{
	cerr<<"EXEC Branch FMD digits not found"<<endl;
	return;
      } 
	  
      if(plFMD->TreeR()==0) plFMD->MakeTree("R");

      //Make branches
      const Int_t kBufferSize = 16000;
      plFMD->TreeR()->Branch("FMD", &reconParticles, kBufferSize);

      
      Int_t zeroADC=6;
      // read Digits 
      AliFMDdigit  *fmdDigit;
       if (plFMD->TreeD()->GetEvent(0))
	{	  	  
	  Int_t nDigits=digits->GetEntries();
	  Int_t recParticles[6];
	  Int_t nRecPart=0 ;
	  Int_t zeroPads=0;
	  Int_t numberOfPads=0;
	  Int_t pedestal;
	  Float_t channelWidth=(22400*50)/1024;
	  for (Int_t digit=0;digit<nDigits;digit++) 
	    {
	      fmdDigit=(AliFMDdigit*)digits->UncheckedAt(digit);    
	      ivol=fmdDigit->Volume();
	      iSector=fmdDigit->NumberOfSector();
	      iRing=fmdDigit->NumberOfRing();
	      pedestal=Int_t(gRandom->Gaus(500,250));
	      padADC= fmdDigit->ADCsignal()
		-Int_t(Float_t(pedestal)/channelWidth);
	      if (padADC<0) padADC=0;
              hTotal[ivol]->Fill(iSector,iRing,padADC);
	    } //digit loop

	  //reconstruct multiplicity in 0.1 eta according Poisson distribution 

	  Int_t rmin=0; Int_t rmax=0; 
	  Int_t smin=0; Int_t smax=0; 
	  for (ivol=0; ivol<knumVolumes; ivol++)
	    {
	      Float_t ring2number=Float_t (numberOfRings[ivol])/
		(rout[ivol]-rin[ivol]);
	      Float_t realZ=zVertex+z[ivol];
	      theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
	      etain = - TMath::Log( TMath::Tan(theta/2.));
	      theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
	      etaout=- TMath::Log( TMath::Tan(theta/2.));
	      numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
	      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));
		  rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
		  eta=eta+0.1;
		  theta = 2.*TMath::ATan (TMath::Exp (-eta));
		  rad = TMath::Abs(realZ) * (TMath::Tan (theta));
		  rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
		  
		  zeroPads=0;
		  smin=0;
		  smax=numberOfSectors[ivol]; 
		  numberOfPads=(rmax-rmin)*(smax-smin);
		  for (Int_t iring=rmin; iring<rmax; iring++) 
		    {
		      for 
			(Int_t isector=0;
			 isector<numberOfSectors[ivol]; 
			 isector++) 
			{
			  if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
			     <zeroADC) zeroPads++;}
		      
		    } //ring
		  Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);

		  Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
		  Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
		  recParticles[0]=ivol+1;
		  recParticles[1]=smin;
		  recParticles[2]=smax;
		  recParticles[3]=rmin;
		  recParticles[4]=rmax;
		  recParticles[5]= fRecon;
		  new((*reconParticles)[nRecPart++])   AliFMDReconstParticles(recParticles);             
		  

		 } // eta
	     } // volume
	   
	}//if (plFMD->TreeD()->GetEvent(0))
       plFMD->TreeR()->Reset();
       plFMD->TreeR()->Fill(); 
       plFMD->WriteRecPoints("OVERWRITE");
       plFMD->UnloadDigits();
    } //event loop
  plFMD->UnloadRecPoints();
#ifdef DEBUG
  Info(" Exec"," finished");
#endif
  //  delete hTotal[10]; 

}


//_____________________________________________________________________________
void AliFMDReconstructor::FillESD(AliRunLoader* /*runLoader*/, 
				  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
	16	//_________________________________________________________________________
121a60bd	17	// This is a class that constructs ReconstParticles (reconstructed particles)
37c55dc0	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"
37c55dc0	29	#include "TFolder.h"
46501dfb	30	#include "TH2F.h"
44cc5e44	31	#include "TRandom.h"
37c55dc0	32
dc8af42e	33	// --- Standard library ---
37c55dc0	34	#include <stdlib.h>
88cb7938	35	#include <Riostream.h>
dc8af42e	36
	37	// --- AliRoot header files ---
	38
88cb7938	39	#include "AliRunLoader.h"
	40	#include "AliLoader.h"
	41
dc8af42e	42	#include "AliFMDdigit.h"
dc8af42e	43	#include "AliFMDReconstParticles.h"
121a60bd	44	#include "AliFMDReconstructor.h"
dc8af42e	45	#include "AliRun.h"
88cb7938	46	#include "AliConfig.h"
1a42d4f2	47	#include "AliHeader.h"
1a42d4f2	48	#include "AliGenEventHeader.h"
88cb7938	49
121a60bd	50	ClassImp(AliFMDReconstructor)
dc8af42e	51
dc8af42e	52
88cb7938	53	//____________________________________________________________________________
dc8af42e	54
121a60bd	55	void AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e	56	{
88cb7938	57	//Collects all digits in the same active volume into number of particles
a39f3926	58	/*
88cb7938	59	Reconstruct number of particles
a39f3926	60	in given group of pads for given FMDvolume
	61	determine by numberOfVolume ,
	62	numberOfMinSector,numberOfMaxSector,
	63	numberOfMinRing, numberOgMaxRing
121a60bd	64	Reconstructor method choose dependence on number of empty pads
a39f3926	65	*/
	66
	67
3d44ce66	68	#ifdef DEBUG
	69	Info("Exec ","Start");
	70	#endif
	71
a39f3926	72	Int_t const knumVolumes=5;
46501dfb	73	Int_t padADC;
1a42d4f2	74	Float_t eta, etain,etaout,rad,theta;
cb1df35e	75	Int_t ivol, iSector, iRing;
1a42d4f2	76	Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
1a42d4f2	77	Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
2ad713c4	78	Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
46501dfb	79	Int_t numberOfRings[5]={512,256,512,256,512};
a39f3926	80	Int_t numberOfSectors[5]= {20,40,20,40,20};
a39f3926	81	Int_t numberOfEtaIntervals[5];
cb1df35e	82	// number of ring for boundary 0.1 eta
cb1df35e	83
88cb7938	84
121a60bd	85	if (runLoader == 0x0)
88cb7938	86	{
	87	Error("Exec","Run Loader loader is NULL - Session not opened");
	88	return;
	89	}
3d44ce66	90
121a60bd	91	AliLoader* plFMD = runLoader->GetLoader("FMDLoader");
3d44ce66	92	if (plFMD == 0x0)
88cb7938	93	{
121a60bd	94	Fatal("AliFMDReconstructor","Can not find FMD (loader) in specified event");
88cb7938	95	return;//never reached
	96	}
	97
121a60bd	98	if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
121a60bd	99	if (!runLoader->TreeE()) runLoader->LoadHeader();
3d44ce66	100
46501dfb	101	TDirectory* cwd = gDirectory;
46501dfb	102	gDirectory = 0x0;
3d44ce66	103	Text_t buf1[20];
46501dfb	104	TH2F* hTotal[10];
	105	for (Int_t j=1; j<=5; j++){
	106	sprintf(buf1,"hTotal%d",j);
	107
	108	hTotal[j] = new TH2F(buf1," Number of primary particles ",
	109	numberOfSectors[j-1],1,numberOfSectors[j-1],
	110	numberOfRings[j-1],1,numberOfRings[j-1]);
	111	}
	112	gDirectory = cwd;
	113
121a60bd	114	plFMD->LoadRecPoints("RECREATE");
44cc5e44	115	TClonesArray* reconParticles = new TClonesArray("AliFMDReconstParticles");
fcc1c32e	116	TClonesArray* digits = new TClonesArray("AliFMDdigit");
44cc5e44	117
3d44ce66	118	Int_t retval=0;
121a60bd	119	Int_t nevents=Int_t (runLoader->TreeE()->GetEntries());
3d44ce66	120	#ifdef DEBUG
121a60bd	121	cout<<" nevents "<<nevents<<endl;
3d44ce66	122	#endif
121a60bd	123	for(Int_t ievent=0;ievent<nevents;ievent++)
dc8af42e	124	{
3d44ce66	125	#ifdef DEBUG
	126	cout<<" *** event "<<ievent<<endl;
	127	#endif
121a60bd	128	runLoader->GetEvent(ievent) ;
3d44ce66	129	//event z-vertex for correction eta-rad dependence
121a60bd	130	AliHeader *header = runLoader->GetHeader();
3d44ce66	131	AliGenEventHeader* genHeader = header->GenEventHeader();
3d44ce66	132	TArrayF *o = new TArrayF(3);
27301b0a	133	if (genHeader) genHeader->PrimaryVertex(*o);
3d44ce66	134	Float_t zVertex=o->At(2);
88cb7938	135
3d44ce66	136	for (Int_t i=0; i<5; i++)
46501dfb	137	hTotal[i+1]->Reset();
46501dfb	138
3d44ce66	139	retval = plFMD->LoadDigits("READ");
46501dfb	140	if (retval)
	141	{
	142	Error("Exec","Error occured while loading digits. Exiting.");
	143	return;
	144	}
	145
3d44ce66	146	TTree* treeD = plFMD->TreeD();
46501dfb	147	if (treeD == 0x0)
	148	{
	149	Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
	150	return;
	151	}
	152
	153	TBranch *brDigits=0;
	154
	155	brDigits=treeD->GetBranch("FMD");
3d44ce66	156
46501dfb	157	if (brDigits) {
44cc5e44	158	brDigits->SetAddress(&digits);
46501dfb	159	}else{
	160	cerr<<"EXEC Branch FMD digits not found"<<endl;
	161	return;
	162	}
	163
3d44ce66	164	if(plFMD->TreeR()==0) plFMD->MakeTree("R");
88cb7938	165
dc8af42e	166	//Make branches
44cc5e44	167	const Int_t kBufferSize = 16000;
44cc5e44	168	plFMD->TreeR()->Branch("FMD", &reconParticles, kBufferSize);
1a42d4f2	169
37c55dc0	170
46501dfb	171	Int_t zeroADC=6;
3d44ce66	172	// read Digits
1a42d4f2	173	AliFMDdigit *fmdDigit;
44cc5e44	174	if (plFMD->TreeD()->GetEvent(0))
	175	{
	176	Int_t nDigits=digits->GetEntries();
46501dfb	177	Int_t recParticles[6];
46501dfb	178	Int_t nRecPart=0 ;
3d44ce66	179	Int_t zeroPads=0;
	180	Int_t numberOfPads=0;
	181	Int_t pedestal;
	182	Float_t channelWidth=(22400*50)/1024;
	183	for (Int_t digit=0;digit<nDigits;digit++)
	184	{
44cc5e44	185	fmdDigit=(AliFMDdigit*)digits->UncheckedAt(digit);
3d44ce66	186	ivol=fmdDigit->Volume();
	187	iSector=fmdDigit->NumberOfSector();
	188	iRing=fmdDigit->NumberOfRing();
	189	pedestal=Int_t(gRandom->Gaus(500,250));
	190	padADC= fmdDigit->ADCsignal()
	191	-Int_t(Float_t(pedestal)/channelWidth);
	192	if (padADC<0) padADC=0;
46501dfb	193	hTotal[ivol]->Fill(iSector,iRing,padADC);
3d44ce66	194	} //digit loop
	195
	196	//reconstruct multiplicity in 0.1 eta according Poisson distribution
46501dfb	197
3d44ce66	198	Int_t rmin=0; Int_t rmax=0;
	199	Int_t smin=0; Int_t smax=0;
	200	for (ivol=0; ivol<knumVolumes; ivol++)
	201	{
	202	Float_t ring2number=Float_t (numberOfRings[ivol])/
	203	(rout[ivol]-rin[ivol]);
	204	Float_t realZ=zVertex+z[ivol];
	205	theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
	206	etain = - TMath::Log( TMath::Tan(theta/2.));
	207	theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
	208	etaout=- TMath::Log( TMath::Tan(theta/2.));
	209	numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
	210	eta=etain;
	211	for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
	212	{
	213	theta = 2.*TMath::ATan (TMath::Exp (-eta));
	214	Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
	215	rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
	216	eta=eta+0.1;
	217	theta = 2.*TMath::ATan (TMath::Exp (-eta));
	218	rad = TMath::Abs(realZ) * (TMath::Tan (theta));
	219	rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
	220
	221	zeroPads=0;
	222	smin=0;
	223	smax=numberOfSectors[ivol];
	224	numberOfPads=(rmax-rmin)*(smax-smin);
	225	for (Int_t iring=rmin; iring<rmax; iring++)
	226	{
	227	for
	228	(Int_t isector=0;
	229	isector<numberOfSectors[ivol];
	230	isector++)
	231	{
	232	if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
	233	<zeroADC) zeroPads++;}
	234
	235	} //ring
	236	Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
	237
	238	Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
	239	Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
	240	recParticles[0]=ivol+1;
	241	recParticles[1]=smin;
	242	recParticles[2]=smax;
	243	recParticles[3]=rmin;
	244	recParticles[4]=rmax;
	245	recParticles[5]= fRecon;
44cc5e44	246	new((*reconParticles)[nRecPart++]) AliFMDReconstParticles(recParticles);
3d44ce66	247
46501dfb	248
1a42d4f2	249	} // eta
1a42d4f2	250	} // volume
cb1df35e	251
44cc5e44	252	}//if (plFMD->TreeD()->GetEvent(0))
3d44ce66	253	plFMD->TreeR()->Reset();
	254	plFMD->TreeR()->Fill();
	255	plFMD->WriteRecPoints("OVERWRITE");
	256	plFMD->UnloadDigits();
dc8af42e	257	} //event loop
121a60bd	258	plFMD->UnloadRecPoints();
3d44ce66	259	#ifdef DEBUG
	260	Info(" Exec"," finished");
	261	#endif
	262	// delete hTotal[10];
	263
88cb7938	264	}
dc8af42e	265
121a60bd	266
	267	//_____________________________________________________________________________
	268	void AliFMDReconstructor::FillESD(AliRunLoader* /runLoader/,
	269	AliESD* /esd/) const
	270	{
	271	// nothing to be done
	272
	273	}
	274