[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>

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