[u/mrichter/AliRoot.git] / FIT / AliFITReconstructor.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.                  *
 **************************************************************************/

/*********************************************************************
 *  FIT reconstruction and filling ESD
 *  - reconstruct mean time (interation time) 
 *  - vertex position
 *  -  multiplicity
 ********************************************************************/

#include <AliESDEvent.h>
#include "AliLog.h"
#include "AliFITRecPoint.h"
#include "AliFITDigit.h"
#include "AliFITReconstructor.h"
#include "AliESDFIT.h"
#include "AliLog.h"
#include "AliESDRun.h"
#include "AliFITRawReader.h"

#include <TArrayI.h>
#include <TGraph.h>
#include <TMath.h>
#include <Riostream.h>

using std::cout;
using std::endl;

ClassImp(AliFITReconstructor)

AliFITReconstructor:: AliFITReconstructor(): AliReconstructor(),
					     fESD(NULL),
					     fESDFIT(NULL),
					     fDigits(NULL)

{
 
  //constructor
 
  
   fESDFIT  = new AliESDFIT();

 
}
//_____________________________________________________________________________
void AliFITReconstructor::Init()
{
// initializer

  fESDFIT  = new AliESDFIT;
 }

//______________________________________________________________________
void AliFITReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
{
// converts RAW to digits 
  cout<<" AliFITReconstructor::ConvertDigits "<<rawReader<<" "<<digitsTree<<endl;
  if (!digitsTree) {
    AliError("No digits tree!");
    return;
  }

  fDigits = new TClonesArray ("AliFITDigit", 100);
  digitsTree->Branch("FIT", &fDigits);

  AliFITRawReader myrawreader(rawReader);
  if (myrawreader.Next()) 
    {
      Int_t allData[500];
      for (Int_t i=0; i<500; i++)  allData[i]=0;
      for (Int_t i=0; i<500; i++) 
	if(myrawreader.GetData(i)>0)  { 
	  allData[i]=myrawreader.GetData(i);
	}
   
      Int_t timeCFD, timeLED, timeQT1, timeQT0;
      for (Int_t ipmt=0; ipmt<160; ipmt++) {
	if(allData[ipmt]>0) {
	  timeCFD = allData[ipmt];
	  timeLED = allData[ipmt];
	  timeQT0= allData[ipmt+160];
	  timeQT1 = allData[ipmt+320];
	  //add digit
	  new((*fDigits)[fDigits->GetEntriesFast()] )
	    AliFITDigit( ipmt,timeCFD, timeLED, timeQT0, timeQT1, 0);
	}
      }
      
      
      
      digitsTree->Fill(); 
    }

  
}

 //____________________________________________________________
void AliFITReconstructor::FillESD(TTree *digitsTree, TTree * /*clustersTree*/, AliESDEvent *pESD) const
{
  
  /***************************************************
  Resonstruct digits to vertex position
  ****************************************************/
  
  AliDebug(1,Form("Start FillESD FIT"));
  if(!pESD) {
    AliError("No ESD Event");
    return;
  }
  
  Float_t channelWidth = 24.4;  
  Float_t c = 0.0299792458; // cm/ps
  Float_t currentVertex, shift=0;
  Int_t ncont=-1;
   const AliESDVertex* vertex = pESD->GetPrimaryVertex();
  if (!vertex)        vertex = pESD->GetPrimaryVertexSPD();
  if (!vertex)        vertex = pESD->GetPrimaryVertexTPC();
  if (!vertex)        vertex = pESD->GetVertex();
  
  if (vertex) {
    AliDebug(2, Form("Got %s (%s) from ESD: %f", 
		     vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
    currentVertex = vertex->GetZ();
    
    ncont = vertex->GetNContributors();
    if(ncont>0 ) {
      shift = currentVertex/c;
    }
  } //vertex

 
  // FIT digits reconstruction
  
  if (!digitsTree) {
    AliError("No digits tree!");
    return;
  }
  
  fDigits = new TClonesArray("AliFITDigit", 100);
  // digitsTree->Print();
  TBranch *brDigits=digitsTree->GetBranch("FIT");
  if (brDigits) 
    brDigits->SetAddress(&fDigits);
  else {
    AliError("no FIT branch");
    return;
  }
  const Float_t max_time = 1e7;
  Int_t pmt=-1;
  Float_t  time[160],amp[160];  
  for(Int_t i=0; i<160; i++)   {  time[i]=max_time; amp[i]=0;}
  
  Int_t nEntries = (Int_t)digitsTree->GetEntries();
  for (Int_t iev=0; iev<nEntries; iev++) {
    digitsTree->GetEvent(iev,1);
    brDigits->GetEvent(iev);
    Int_t nDigits = fDigits->GetEntriesFast(); 
    for (Int_t dig=0; dig<nDigits; dig++) {    
      AliFITDigit* digit = (AliFITDigit*) fDigits->At(dig);      
      pmt = digit->NPMT();
      time[pmt] = Float_t (digit->TimeCFD() );
      amp[pmt] = Float_t (digit->TimeQT1() - digit->TimeQT0() );
    } 
    fESDFIT->SetFITtime(time);         // best TOF on each PMT 
    fESDFIT->SetFITamplitude(amp);     // number of particles(MIPs) on each 
    
    Float_t firsttime[3] = {max_time,max_time,max_time};
    
    Float_t vertexFIT = 9999999;
    for (Int_t ipmt=0; ipmt<80; ipmt++)//timeC
 	  if(time[ipmt]<firsttime[2]) firsttime[2]=time[ipmt]; 
    
    for ( Int_t ipmt=80; ipmt<160; ipmt++) 
	  if(time[ipmt]<firsttime[1]) firsttime[1]=time[ipmt]; 
    
    
    AliDebug(5,Form("1stimeA %f , 1sttimeC %f  ",
		    firsttime[1],
		    firsttime[2] ) );
    if (firsttime[1]<max_time && firsttime[2]<max_time)  {
	firsttime[0] =  channelWidth *(firsttime[1] + firsttime[2])/2;
	vertexFIT =  c*channelWidth*(firsttime[1] - firsttime[2])/2;
    }
      if (firsttime[1]<max_time) 
	firsttime[1] = firsttime[1] * channelWidth + shift; 
      
      if (firsttime[2]<max_time) 
	firsttime[2] = firsttime[2] * channelWidth - shift; 
      
      
      
      for(Int_t i=0; i<3; i++) {
	fESDFIT->SetFITT0(i,firsttime[i]);   // interaction time (ns) 
      }
      fESDFIT->SetFITzVertex(vertexFIT); //vertex Z position 
      
      
      AliDebug(1,Form("FIT: SPDshift %f Vertex %f  FITsignal %f ps FITA %f ps FITC %f ps \n",shift, vertexFIT, firsttime[0], firsttime[1],firsttime[2]));
  }
  
   }
Commit	Line	Data
c1c44db3	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	* FIT reconstruction and filling ESD
	18	* - reconstruct mean time (interation time)
	19	* - vertex position
	20	* - multiplicity
	21	********************************************************************/
	22
	23	#include <AliESDEvent.h>
	24	#include "AliLog.h"
	25	#include "AliFITRecPoint.h"
	26	#include "AliFITDigit.h"
	27	#include "AliFITReconstructor.h"
	28	#include "AliESDFIT.h"
	29	#include "AliLog.h"
	30	#include "AliESDRun.h"
	31	#include "AliFITRawReader.h"
	32
	33	#include <TArrayI.h>
	34	#include <TGraph.h>
	35	#include <TMath.h>
	36	#include <Riostream.h>
	37
5e2539d9	38	using std::cout;
	39	using std::endl;
	40
c1c44db3	41	ClassImp(AliFITReconstructor)
	42
	43	AliFITReconstructor:: AliFITReconstructor(): AliReconstructor(),
	44	fESD(NULL),
292cfeb6	45	fESDFIT(NULL),
292cfeb6	46	fDigits(NULL)
c1c44db3	47
	48	{
	49
	50	//constructor
	51
	52
	53	fESDFIT = new AliESDFIT();
	54
	55
	56	}
	57	//_____________________________________________________________________________
	58	void AliFITReconstructor::Init()
	59	{
	60	// initializer
	61
	62	fESDFIT = new AliESDFIT;
	63	}
	64
	65	//______________________________________________________________________
	66	void AliFITReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
	67	{
	68	// converts RAW to digits
	69	cout<<" AliFITReconstructor::ConvertDigits "<<rawReader<<" "<<digitsTree<<endl;
	70	if (!digitsTree) {
	71	AliError("No digits tree!");
	72	return;
	73	}
	74
	75	fDigits = new TClonesArray ("AliFITDigit", 100);
	76	digitsTree->Branch("FIT", &fDigits);
	77
	78	AliFITRawReader myrawreader(rawReader);
	79	if (myrawreader.Next())
	80	{
	81	Int_t allData[500];
	82	for (Int_t i=0; i<500; i++) allData[i]=0;
	83	for (Int_t i=0; i<500; i++)
	84	if(myrawreader.GetData(i)>0) {
	85	allData[i]=myrawreader.GetData(i);
	86	}
	87
	88	Int_t timeCFD, timeLED, timeQT1, timeQT0;
	89	for (Int_t ipmt=0; ipmt<160; ipmt++) {
	90	if(allData[ipmt]>0) {
	91	timeCFD = allData[ipmt];
	92	timeLED = allData[ipmt];
	93	timeQT0= allData[ipmt+160];
	94	timeQT1 = allData[ipmt+320];
	95	//add digit
	96	new((*fDigits)[fDigits->GetEntriesFast()] )
	97	AliFITDigit( ipmt,timeCFD, timeLED, timeQT0, timeQT1, 0);
	98	}
	99	}
	100
	101
	102
	103	digitsTree->Fill();
	104	}
	105
	106
	107	}
	108
	109	//____________________________________________________________
292cfeb6	110	void AliFITReconstructor::FillESD(TTree digitsTree, TTree /clustersTree/, AliESDEvent *pESD) const
c1c44db3	111	{
	112
	113	/***************************************************
	114	Resonstruct digits to vertex position
	115	****************************************************/
	116
	117	AliDebug(1,Form("Start FillESD FIT"));
	118	if(!pESD) {
	119	AliError("No ESD Event");
	120	return;
	121	}
	122
	123	Float_t channelWidth = 24.4;
	124	Float_t c = 0.0299792458; // cm/ps
292cfeb6	125	Float_t currentVertex, shift=0;
c1c44db3	126	Int_t ncont=-1;
292cfeb6	127	const AliESDVertex* vertex = pESD->GetPrimaryVertex();
c1c44db3	128	if (!vertex) vertex = pESD->GetPrimaryVertexSPD();
	129	if (!vertex) vertex = pESD->GetPrimaryVertexTPC();
	130	if (!vertex) vertex = pESD->GetVertex();
	131
	132	if (vertex) {
	133	AliDebug(2, Form("Got %s (%s) from ESD: %f",
	134	vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
	135	currentVertex = vertex->GetZ();
	136
	137	ncont = vertex->GetNContributors();
	138	if(ncont>0 ) {
	139	shift = currentVertex/c;
	140	}
	141	} //vertex
	142
292cfeb6	143
c1c44db3	144	// FIT digits reconstruction
	145
	146	if (!digitsTree) {
	147	AliError("No digits tree!");
	148	return;
	149	}
	150
	151	fDigits = new TClonesArray("AliFITDigit", 100);
	152	// digitsTree->Print();
	153	TBranch *brDigits=digitsTree->GetBranch("FIT");
	154	if (brDigits)
	155	brDigits->SetAddress(&fDigits);
	156	else {
	157	AliError("no FIT branch");
	158	return;
	159	}
	160	const Float_t max_time = 1e7;
	161	Int_t pmt=-1;
	162	Float_t time[160],amp[160];
	163	for(Int_t i=0; i<160; i++) { time[i]=max_time; amp[i]=0;}
	164
	165	Int_t nEntries = (Int_t)digitsTree->GetEntries();
	166	for (Int_t iev=0; iev<nEntries; iev++) {
	167	digitsTree->GetEvent(iev,1);
	168	brDigits->GetEvent(iev);
	169	Int_t nDigits = fDigits->GetEntriesFast();
	170	for (Int_t dig=0; dig<nDigits; dig++) {
	171	AliFITDigit* digit = (AliFITDigit*) fDigits->At(dig);
	172	pmt = digit->NPMT();
	173	time[pmt] = Float_t (digit->TimeCFD() );
	174	amp[pmt] = Float_t (digit->TimeQT1() - digit->TimeQT0() );
	175	}
	176	fESDFIT->SetFITtime(time); // best TOF on each PMT
	177	fESDFIT->SetFITamplitude(amp); // number of particles(MIPs) on each
	178
	179	Float_t firsttime[3] = {max_time,max_time,max_time};
	180
	181	Float_t vertexFIT = 9999999;
	182	for (Int_t ipmt=0; ipmt<80; ipmt++)//timeC
	183	if(time[ipmt]<firsttime[2]) firsttime[2]=time[ipmt];
	184
	185	for ( Int_t ipmt=80; ipmt<160; ipmt++)
	186	if(time[ipmt]<firsttime[1]) firsttime[1]=time[ipmt];
	187
	188
	189	AliDebug(5,Form("1stimeA %f , 1sttimeC %f ",
	190	firsttime[1],
	191	firsttime[2] ) );
	192	if (firsttime[1]<max_time && firsttime[2]<max_time) {
	193	firsttime[0] = channelWidth *(firsttime[1] + firsttime[2])/2;
	194	vertexFIT = cchannelWidth(firsttime[1] - firsttime[2])/2;
	195	}
	196	if (firsttime[1]<max_time)
	197	firsttime[1] = firsttime[1] * channelWidth + shift;
	198
	199	if (firsttime[2]<max_time)
	200	firsttime[2] = firsttime[2] * channelWidth - shift;
	201
	202
	203
	204	for(Int_t i=0; i<3; i++) {
	205	fESDFIT->SetFITT0(i,firsttime[i]); // interaction time (ns)
	206	}
	207	fESDFIT->SetFITzVertex(vertexFIT); //vertex Z position
208
209
210	AliDebug(1,Form("FIT: SPDshift %f Vertex %f FITsignal %f ps FITA %f ps FITC %f ps \n",shift, vertexFIT, firsttime[0], firsttime[1],firsttime[2]));
211	}
212
213	}