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

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

#include <AliESDEvent.h>
#include "AliLog.h"
#include "AliT0RecPoint.h"
#include "AliRawReader.h"
#include "AliT0RawReader.h"
#include "AliT0digit.h"
#include "AliT0Reconstructor.h"
#include "AliT0Parameters.h"
#include "AliT0Calibrator.h"
#include "AliESDfriend.h"
#include "AliESDTZEROfriend.h"
#include "AliLog.h"

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

ClassImp(AliT0Reconstructor)

  AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
					     fdZonA(0),
					     fdZonC(0),
					     fZposition(0),
					     fParam(NULL),
					     fAmpLEDrec(),
					     fQTC(0),
					     fAmpLED(0),
                                             fCalib(),
                                             fLatencyHPTDC(9000),
                                             fLatencyL1(0),
                                             fLatencyL1A(0),
                                             fLatencyL1C(0),
                                             fESDTZEROfriend(NULL)

{
  //constructor

  fParam = AliT0Parameters::Instance();
  fParam->Init();
 
  for (Int_t i=0; i<24; i++){
        TGraph* gr = fParam ->GetAmpLEDRec(i);
	if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ; 
	  TGraph* gr1 = fParam ->GetAmpLED(i);
	  if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ; 
	  TGraph* gr2 = fParam ->GetQTC(i);
	  if (gr2) fQTC.AddAtAndExpand(gr2,i) ; 	
  }

  fLatencyL1 = fParam->GetLatencyL1();
  fLatencyL1A = fParam->GetLatencyL1A();
  fLatencyL1C = fParam->GetLatencyL1C();
  fLatencyHPTDC = fParam->GetLatencyHPTDC();
  AliDebug(10,Form(" LatencyL1 %f latencyL1A %f latencyL1C %f latencyHPTDC %f \n",fLatencyL1, fLatencyL1A, fLatencyL1C, fLatencyHPTDC));
  
  // fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
  //fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
  //here real Z position
  fdZonC = TMath::Abs(fParam->GetZPosition("T0/C/PMT1"));
  fdZonA = TMath::Abs(fParam->GetZPosition("T0/A/PMT15"));

  fCalib = new AliT0Calibrator();
  fESDTZEROfriend = new AliESDTZEROfriend();

}

//_____________________________________________________________________________
void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
  
{
  // T0 digits reconstruction
  Int_t refAmp = Int_t (GetRecoParam()->GetRefAmp());
  
  TArrayI * timeCFD = new TArrayI(24); 
  TArrayI * timeLED = new TArrayI(24); 
  TArrayI * chargeQT0 = new TArrayI(24); 
  TArrayI * chargeQT1 = new TArrayI(24); 

 
  Float_t channelWidth = fParam->GetChannelWidth() ;  
  Float_t meanVertex = fParam->GetMeanVertex();
  Float_t c = 0.0299792; // cm/ps
  Double32_t vertex = 9999999;
  Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;

  
  AliDebug(1,Form("Start DIGITS reconstruction "));
  

  TBranch *brDigits=digitsTree->GetBranch("T0");
  AliT0digit *fDigits = new AliT0digit() ;
  if (brDigits) {
    brDigits->SetAddress(&fDigits);
  }else{
    AliError(Form("EXEC Branch T0 digits not found"));
     return;
  }
  
  digitsTree->GetEvent(0);
  digitsTree->GetEntry(0);
  brDigits->GetEntry(0);
  fDigits->GetTimeCFD(*timeCFD);
  fDigits->GetTimeLED(*timeLED);
  fDigits->GetQT0(*chargeQT0);
  fDigits->GetQT1(*chargeQT1);
  Int_t onlineMean =  fDigits->MeanTime();

  Bool_t tr[5];
  for (Int_t i=0; i<5; i++) tr[i]=false; 
  
  Double32_t besttimeA=999999;
  Double32_t besttimeC=999999;
  Int_t pmtBestA=99999;
  Int_t pmtBestC=99999;
  
  AliT0RecPoint* frecpoints= new AliT0RecPoint ();
  clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
  
  Float_t time[24], adc[24];
  for (Int_t ipmt=0; ipmt<24; ipmt++) {
    if(timeCFD->At(ipmt)>0 ){
     if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)  
	adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
      else
	adc[ipmt] = 0;
      
     time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]),  timeCFD->At(ipmt)) ;
	     
      Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
      //    time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl),  timeCFD->At(ipmt) ) ;
      AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
		       ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));

      Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
      Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
      AliDebug(10,Form("  Amlitude in MIPS LED %f ,  QTC %f in channels %f\n ",ampMip,qtMip, adc[ipmt]));
      
      frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
      frecpoints->SetAmpLED(ipmt, Float_t( ampMip)); //for cosmic &pp beam 
      frecpoints->SetAmp(ipmt, Float_t(qtMip));
      
    }
    else {
      time[ipmt] = 0;
      adc[ipmt] = 0;
    }
  }
  
  for (Int_t ipmt=0; ipmt<12; ipmt++){
    if(time[ipmt] > 1 ) {
      if(time[ipmt]<besttimeC){
	besttimeC=time[ipmt]; //timeC
	pmtBestC=ipmt;
      }
    }
  }
  for ( Int_t ipmt=12; ipmt<24; ipmt++){
    if(time[ipmt] > 1) {
      if(time[ipmt]<besttimeA) {
	besttimeA=time[ipmt]; //timeA
        pmtBestA=ipmt;}
    }
  }
  if(besttimeA < 999999) {
    frecpoints->SetTimeBestA(Int_t(besttimeA *channelWidth - fdZonA/c));
    tr[1]=true;
  }
  if( besttimeC < 999999 ) {
    frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth - fdZonA/c));
    tr[2]=true;
  }
  AliDebug(10,Form(" besttimeA %f ch,  besttimeC %f ch",besttimeA, besttimeC));
  if(besttimeA <999999 && besttimeC < 999999 ){
    //    timeDiff = (besttimeC - besttimeA)*channelWidth;
    timeDiff = (besttimeA - besttimeC)*channelWidth;
    meanTime = (besttimeA + besttimeC)/2;// * channelWidth); 
    timeclock = meanTime *channelWidth -fdZonA/c ;
    vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
    tr[0]=true; 
  }
  frecpoints->SetVertex(vertex);
  frecpoints->SetMeanTime(meanTime);
  frecpoints->SetT0clock(timeclock);
  frecpoints->SetT0Trig(tr);

  AliDebug(10,Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));

  //online mean
  frecpoints->SetOnlineMean(Int_t(onlineMean));
  AliDebug(10,Form("  timeDiff %f #channel,  meanTime %f #channel, vertex %f cm online mean %i timeclock %f ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
  
  
  clustersTree->Fill();

  delete timeCFD;
  delete timeLED;
  delete chargeQT0; 
  delete chargeQT1; 
}


//_______________________________________________________________________

void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
{
  // T0 raw ->
  //
  // reference amplitude and time ref. point from reco param

  Float_t refAmp = GetRecoParam()->GetRefAmp();
  Int_t refPoint = 0;
  //Bad channel
  Int_t badpmt = GetRecoParam()->GetRefPoint();
  
  Int_t low[110], high[110];

  Int_t allData[110][5];
  
  Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
  Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
  Float_t c = 29.9792458; // cm/ns
  Double32_t vertex = 9999999;
  Int_t onlineMean=0;
  // Float_t meanVertex = fParam->GetMeanVertex();
  Float_t meanVertex = 0;
  for (Int_t i0=0; i0<105; i0++)
    {
      for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0; 
      low[i0] = Int_t (GetRecoParam()->GetLow(i0));	
      high[i0] = Int_t (GetRecoParam()->GetHigh(i0));
    }
   
  Double32_t besttimeA=9999999;
  Double32_t besttimeC=9999999;
  Int_t pmtBestA=99999;
  Int_t pmtBestC=99999;
   
  AliT0RecPoint* frecpoints= new AliT0RecPoint ();
  
  recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
   
  AliDebug(10," before read data ");
  AliT0RawReader myrawreader(rawReader);

  UInt_t type =rawReader->GetType();

  if (!myrawreader.Next())
    AliDebug(1,Form(" no raw data found!!"));
  else
    {  
   for (Int_t i=0; i<24; i++)
    {
      timeCFD[i]=0; timeLED[i]=0; chargeQT0[i]=0; chargeQT1[i]=0;
    }
     Int_t fBCID=Int_t (rawReader->GetBCID());
      Int_t trmbunch= myrawreader.GetTRMBunchID();
      if(type == 7  ) {  //only physics 
	for (Int_t i=0; i<105; i++) {
	for (Int_t iHit=0; iHit<5; iHit++) 
	  {
	    allData[i][iHit] = myrawreader.GetData(i,iHit);
	  }
	}
	Int_t ref=0;
	if (refPoint>0) 
	  ref = allData[refPoint][0]-5000;
	
	Float_t channelWidth = fParam->GetChannelWidth() ;  
	
	//       Int_t meanT0 = fParam->GetMeanT0();
	
	    for (Int_t in=0; in<12; in++)  
	      {
		for (Int_t iHit=0; iHit<5; iHit++) 
		  {
		    if(allData[in+1][iHit] > low[in+1] && 
		       allData[in+1][iHit] < high[in+1])
		      {
			timeCFD[in] = allData[in+1][iHit] ; 
			break;
		      }
		  }
		for (Int_t iHit=0; iHit<5; iHit++) 
		  {
		    if(allData[in+12+1][iHit] > low[in+12+1] && 
		       allData[in+1+12][iHit] < high[in+12+1])
		      {
			timeLED[in] = allData[in+12+1][iHit] ;
			break;
		      }
		  }
		for (Int_t iHit=0; iHit<5; iHit++) 
		  {
		    if(allData[in+1+56][iHit] > low[in+1+56] && 
		       allData[in+1+56][iHit] < high[in+1+56])
		      {
			timeCFD[in+12] = allData[in+56+1][iHit] ;
			break;
		      }
		  }
		
		for (Int_t iHit=0; iHit<5; iHit++) 
		  {
		    if(allData[in+1+68][iHit] > low[in+1+68] && 
		       allData[in+1+68][iHit] < high[in+1+68])
		      {
			timeLED[in+12] = allData[in+68+1][iHit] ;
			break;
		      }
		  }
		AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
				  in, timeCFD[in],timeCFD[in+12],timeLED[in], 
				  timeLED[in+12]));   
		
	      }
	    
	  
	for (Int_t iHit=0; iHit<5; iHit++) 
	  {
	    for (Int_t in=0; in<12;  in++)
	      {
		if(allData[2*in+25][iHit] > low[2*in+25] && 
		   allData[2*in+25][iHit] < high[2*in+25])
		  chargeQT0[in]=allData[2*in+25][iHit];
		if(allData[2*in+26][iHit] > low[2*in+26] && 
		   allData[2*in+26][iHit] < high[2*in+26])
		  chargeQT1[in]=allData[2*in+26][iHit];
		
	      }
	    
	    for (Int_t in=12; in<24;  in++)
	      {
		if(allData[2*in+57][iHit] > low[2*in+57] && 
		   allData[2*in+57][iHit] < high[2*in+57])		  
		  chargeQT0[in]=allData[2*in+57][iHit];
		if(allData[2*in+58][iHit] > low[2*in+58] && 
		   allData[2*in+58][iHit] < high[2*in+58])		  
		  chargeQT1[in]=allData[2*in+58][iHit];
		
	      }
	  }
	
	
	for (Int_t in=0; in<24; in++)  
	  AliDebug(10, Form(" readed Raw %i %i %i %i %i",
			    in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
   	for (Int_t iHit=0; iHit<5; iHit++) 
	  {
		if(allData[49][iHit] > low[49] && 
		   allData[49][iHit] < high[49])
		  onlineMean = allData[49][iHit];       
	  }
	Double32_t time[24], adc[24],  noncalibtime[24];
	for (Int_t ipmt=0; ipmt<24; ipmt++) {
	  if(timeCFD[ipmt]>0 ){
	   //for simulated data
	     //for physics  data
	   if(( chargeQT0[ipmt] - chargeQT1[ipmt])>0)  {
	     adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
	   }
	   else
	     adc[ipmt] = 0;
	   time[ipmt] = fCalib-> WalkCorrection(Int_t (refAmp), ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
	   
      	   Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
	   // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
	   AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
			    ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
	   Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	   Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
	   AliDebug(10,Form("  Amlitude in MIPS LED %f ; QTC %f;  in channels %f\n ",ampMip,qtMip, adc[ipmt]));
	   //bad peak removing
	   if(sl<550) {
	     frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	     // frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
	     frecpoints->SetAmp(ipmt, Double32_t( qtMip)); //for cosmic &pp beam 
	     frecpoints->SetAmpLED(ipmt, Double32_t(ampMip));	     
	     noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
	   }
	 }
	 else {
	   time[ipmt] = 0;
	   adc[ipmt] = 0;
	   noncalibtime[ipmt] = 0;
	 }
       }
       fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;     
       for (Int_t ipmt=0; ipmt<12; ipmt++){
	 if(time[ipmt] > 1 && ipmt != badpmt )
	   {
	     if(time[ipmt]<besttimeC){
	       besttimeC=time[ipmt]; //timeC
	       pmtBestC=ipmt;
	     }
	   }
       }
       for ( Int_t ipmt=12; ipmt<24; ipmt++)
	 {
	   if(time[ipmt] > 1 && ipmt != badpmt)
	     {
	       if(time[ipmt]<besttimeA) {
		 besttimeA=time[ipmt]; //timeA
		 pmtBestA=ipmt;
	       }
	     }
	 }
       if(besttimeA < 999999) 
	 frecpoints->SetTimeBestA(besttimeA * channelWidth - 1000.*fLatencyHPTDC + 1000.*fLatencyL1A);
       if( besttimeC < 999999 ) 
	 frecpoints->SetTimeBestC(besttimeC * channelWidth - 1000.*fLatencyHPTDC +1000.*fLatencyL1C);
       AliDebug(10,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f ps",
		       pmtBestA,besttimeA, pmtBestC,  besttimeC));
        if(besttimeA <999999 && besttimeC < 999999 ){
	 timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
	 timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - 1000.*fLatencyHPTDC + 1000.*fLatencyL1;  
	 meanTime = (besttimeA+besttimeC-2.*Float_t(ref))/2.;
	 vertex =  meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2; 
	}
      }  //if phys event       
      AliDebug(5,Form("  timeDiff %f #channel,  meanTime %f #channel, TOFmean%f  vertex %f cm meanVertex %f online mean %i \n",timeDiff, meanTime,timeclock, vertex,meanVertex, onlineMean));
      frecpoints->SetT0clock(timeclock);
      frecpoints->SetVertex(vertex);
      frecpoints->SetMeanTime(meanTime);
      frecpoints->SetOnlineMean(Int_t(onlineMean));
	// Set triggers
      
      Bool_t tr[5];
      Int_t trchan[5]= {50,51,52,55,56};
      for (Int_t i=0; i<5; i++) tr[i]=false; 
      for (Int_t itr=0; itr<5; itr++) {
 	for (Int_t iHit=0; iHit<5; iHit++) 
	  {
	    Int_t trr=trchan[itr];
	    if(allData[trr][iHit] > low[trr] &&
	       allData[trr][iHit] < high[trr]&&
	       allData[trr][iHit] > 0) tr[itr]=true;
	  }
      }
      frecpoints->SetT0Trig(tr);
     } // if (else )raw data
  recTree->Fill();
  if(frecpoints) delete frecpoints;
}
  
  
  //____________________________________________________________
  
  void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
  {

  /***************************************************
  Resonstruct digits to vertex position
  ****************************************************/
  
  AliDebug(1,Form("Start FillESD T0"));
  Float_t channelWidth = fParam->GetChannelWidth() ;  
  Float_t c = 0.0299792458; // cm/ps
  Float_t currentVertex=0, shift=0;
  Int_t ncont=0;
  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();
    // cout<<" spdver "<<spdver<<" ncont "<<ncont<<endl;
    if(ncont>2 ) {
      shift = currentVertex/c;
      //	  cout<<" vertex shif "<<shift<<" vertex "<<spdver<<" IsFromVertexer3D  "<<fverSPD->IsFromVertexer3D()<<endl;
    }
  }
  TTree *treeR = clustersTree;
  
  AliT0RecPoint* frecpoints= new AliT0RecPoint ();
  if (!frecpoints) {
    AliError("Reconstruct Fill ESD >> no recpoints found");
    return;
  }
  
  AliDebug(1,Form("Start FillESD T0"));
  TBranch *brRec = treeR->GetBranch("T0");
  if (brRec) {
    brRec->SetAddress(&frecpoints);
  }else{
    AliError(Form("EXEC Branch T0 rec not found"));
    return;
  } 
  
  brRec->GetEntry(0);
  Double32_t amp[24], time[24], ampQTC[24], timecorr[24];  
  Double32_t* tcorr;
  for(Int_t i=0; i<24; i++) 
    amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;

  Double32_t timeClock[3];
  Double32_t zPosition = frecpoints -> GetVertex();
  Double32_t timeStart = frecpoints -> GetMeanTime();
  timeClock[0] = frecpoints -> GetT0clock() ;
  timeClock[1] = frecpoints -> GetBestTimeA() + shift;
  timeClock[2] = frecpoints -> GetBestTimeC() - shift;
  for ( Int_t i=0; i<24; i++) {
    time[i] =  frecpoints -> GetTime(i); // ps to ns
    if ( time[i] >1) {
      ampQTC[i] = frecpoints -> GetAmp(i);
      amp[i] = frecpoints -> AmpLED(i);
    }
  }
  Int_t trig= frecpoints ->GetT0Trig();
  pESD->SetT0Trig(trig);
  
  pESD->SetT0zVertex(zPosition); //vertex Z position 
  pESD->SetT0(timeStart);        // interaction time 
  for(Int_t i=0; i<3; i++) 
    pESD->SetT0TOF(i,timeClock[i]);   // interaction time (ns) 
  pESD->SetT0time(time);         // best TOF on each PMT 
  pESD->SetT0amplitude(ampQTC);     // number of particles(MIPs) on each PMT
  
  AliDebug(1,Form("T0: Vertex %f (T0A+T0C)/2 %f #channels T0signal %f ns OrA %f ns OrC %f T0trig %i\n",zPosition, timeStart, timeClock[0], timeClock[1], timeClock[2], trig));
  
  if (pESD) {
    
    AliESDfriend *fr = (AliESDfriend*)pESD->FindListObject("AliESDfriend");
    if (fr) {
      AliDebug(1, Form("Writing TZERO friend data to ESD tree"));

      if (ncont>2) {
	tcorr = fESDTZEROfriend->GetT0timeCorr();
	for ( Int_t i=0; i<24; i++) {
          timecorr[i]=tcorr[i];
	  if(i<12 && time[i]>1) timecorr[i] -=  shift/channelWidth;
	  if(i>11 && time[i]>1) timecorr[i] +=  shift/channelWidth;
	}
	fESDTZEROfriend->SetT0timeCorr( timecorr) ;
	fESDTZEROfriend->SetT0ampLEDminCFD(amp);
	fESDTZEROfriend->SetT0ampQTC(ampQTC);
	fr->SetTZEROfriend(fESDTZEROfriend);
      }//
    }
  }
     

} // vertex in 3 sigma
Commit	Line	Data
73df58ab	1
dc7ca31d	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	/* $Id$ */
72e48d95	18	/*********************************************************************
	19	* T0 reconstruction and filling ESD
	20	* - reconstruct mean time (interation time)
	21	* - vertex position
	22	* - multiplicity
	23	********************************************************************/
dc7ca31d	24
af885e0f	25	#include <AliESDEvent.h>
dc7ca31d	26	#include "AliLog.h"
dc7ca31d	27	#include "AliT0RecPoint.h"
	28	#include "AliRawReader.h"
	29	#include "AliT0RawReader.h"
dc7ca31d	30	#include "AliT0digit.h"
	31	#include "AliT0Reconstructor.h"
	32	#include "AliT0Parameters.h"
c41ceaac	33	#include "AliT0Calibrator.h"
58bd3a16	34	#include "AliESDfriend.h"
73df58ab	35	#include "AliESDTZEROfriend.h"
8f8d0732	36	#include "AliLog.h"
dc7ca31d	37
	38	#include <TArrayI.h>
	39	#include <TGraph.h>
aad72f45	40	#include <TMath.h>
b09247a2	41	#include <Riostream.h>
dc7ca31d	42
	43	ClassImp(AliT0Reconstructor)
	44
c41ceaac	45	AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
f16935f7	46	fdZonA(0),
	47	fdZonC(0),
	48	fZposition(0),
	49	fParam(NULL),
2e6a5ee0	50	fAmpLEDrec(),
c883fdf2	51	fQTC(0),
c883fdf2	52	fAmpLED(0),
58bd3a16	53	fCalib(),
	54	fLatencyHPTDC(9000),
	55	fLatencyL1(0),
	56	fLatencyL1A(0),
73df58ab	57	fLatencyL1C(0),
73df58ab	58	fESDTZEROfriend(NULL)
58bd3a16	59
e0bba6cc	60	{
72e48d95	61	//constructor
72e48d95	62
74adb36a	63	fParam = AliT0Parameters::Instance();
74adb36a	64	fParam->Init();
c883fdf2	65
74adb36a	66	for (Int_t i=0; i<24; i++){
2e6a5ee0	67	TGraph* gr = fParam ->GetAmpLEDRec(i);
29ed1d0e	68	if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ;
c883fdf2	69	TGraph* gr1 = fParam ->GetAmpLED(i);
	70	if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ;
	71	TGraph* gr2 = fParam ->GetQTC(i);
539b9cb9	72	if (gr2) fQTC.AddAtAndExpand(gr2,i) ;
c883fdf2	73	}
539b9cb9	74
58bd3a16	75	fLatencyL1 = fParam->GetLatencyL1();
	76	fLatencyL1A = fParam->GetLatencyL1A();
	77	fLatencyL1C = fParam->GetLatencyL1C();
	78	fLatencyHPTDC = fParam->GetLatencyHPTDC();
	79	AliDebug(10,Form(" LatencyL1 %f latencyL1A %f latencyL1C %f latencyHPTDC %f \n",fLatencyL1, fLatencyL1A, fLatencyL1C, fLatencyHPTDC));
29ed1d0e	80
adf36b9d	81	// fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
adf36b9d	82	//fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
8f620945	83	//here real Z position
	84	fdZonC = TMath::Abs(fParam->GetZPosition("T0/C/PMT1"));
	85	fdZonA = TMath::Abs(fParam->GetZPosition("T0/A/PMT15"));
539b9cb9	86
12e9daf9	87	fCalib = new AliT0Calibrator();
73df58ab	88	fESDTZEROfriend = new AliESDTZEROfriend();
12e9daf9	89
dc7ca31d	90	}
c41ceaac	91
c41ceaac	92	//_____________________________________________________________________________
dc7ca31d	93	void AliT0Reconstructor::Reconstruct(TTreedigitsTree, TTreeclustersTree) const
94c27e4f	94
dc7ca31d	95	{
94c27e4f	96	// T0 digits reconstruction
38cbfa7c	97	Int_t refAmp = Int_t (GetRecoParam()->GetRefAmp());
776de217	98
c41ceaac	99	TArrayI * timeCFD = new TArrayI(24);
	100	TArrayI * timeLED = new TArrayI(24);
	101	TArrayI * chargeQT0 = new TArrayI(24);
	102	TArrayI * chargeQT1 = new TArrayI(24);
74adb36a	103
d0bcd1fb	104
8955c6b4	105	Float_t channelWidth = fParam->GetChannelWidth() ;
b95e8d87	106	Float_t meanVertex = fParam->GetMeanVertex();
776de217	107	Float_t c = 0.0299792; // cm/ps
adf36b9d	108	Double32_t vertex = 9999999;
adf36b9d	109	Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
776de217	110
94c27e4f	111
dc7ca31d	112	AliDebug(1,Form("Start DIGITS reconstruction "));
94c27e4f	113
2e6a5ee0	114
d0bcd1fb	115	TBranch *brDigits=digitsTree->GetBranch("T0");
e0bba6cc	116	AliT0digit *fDigits = new AliT0digit() ;
dc7ca31d	117	if (brDigits) {
	118	brDigits->SetAddress(&fDigits);
	119	}else{
f16935f7	120	AliError(Form("EXEC Branch T0 digits not found"));
f16935f7	121	return;
dc7ca31d	122	}
e0bba6cc	123
c41ceaac	124	digitsTree->GetEvent(0);
	125	digitsTree->GetEntry(0);
	126	brDigits->GetEntry(0);
	127	fDigits->GetTimeCFD(*timeCFD);
	128	fDigits->GetTimeLED(*timeLED);
	129	fDigits->GetQT0(*chargeQT0);
	130	fDigits->GetQT1(*chargeQT1);
446d6ec4	131	Int_t onlineMean = fDigits->MeanTime();
c883fdf2	132
adf36b9d	133	Bool_t tr[5];
adf36b9d	134	for (Int_t i=0; i<5; i++) tr[i]=false;
c41ceaac	135
adf36b9d	136	Double32_t besttimeA=999999;
adf36b9d	137	Double32_t besttimeC=999999;
c41ceaac	138	Int_t pmtBestA=99999;
c41ceaac	139	Int_t pmtBestC=99999;
dc7ca31d	140
94c27e4f	141	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	142	clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
	143
b95e8d87	144	Float_t time[24], adc[24];
dc7ca31d	145	for (Int_t ipmt=0; ipmt<24; ipmt++) {
c41ceaac	146	if(timeCFD->At(ipmt)>0 ){
d0bcd1fb	147	if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)
d0bcd1fb	148	adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
a7027400	149	else
d0bcd1fb	150	adc[ipmt] = 0;
d0bcd1fb	151
38cbfa7c	152	time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]), timeCFD->At(ipmt)) ;
d0bcd1fb	153
d0bcd1fb	154	Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
8f620945	155	// time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD->At(ipmt) ) ;
d0bcd1fb	156	AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb	157	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
776de217	158
d0bcd1fb	159	Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
d0bcd1fb	160	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
b5a9f753	161	AliDebug(10,Form(" Amlitude in MIPS LED %f , QTC %f in channels %f\n ",ampMip,qtMip, adc[ipmt]));
d0bcd1fb	162
d0bcd1fb	163	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
345f03db	164	frecpoints->SetAmpLED(ipmt, Float_t( ampMip)); //for cosmic &pp beam
345f03db	165	frecpoints->SetAmp(ipmt, Float_t(qtMip));
d0bcd1fb	166
dc7ca31d	167	}
	168	else {
	169	time[ipmt] = 0;
	170	adc[ipmt] = 0;
	171	}
	172	}
94c27e4f	173
dc7ca31d	174	for (Int_t ipmt=0; ipmt<12; ipmt++){
dc7ca31d	175	if(time[ipmt] > 1 ) {
c41ceaac	176	if(time[ipmt]<besttimeC){
	177	besttimeC=time[ipmt]; //timeC
	178	pmtBestC=ipmt;
dc7ca31d	179	}
	180	}
	181	}
	182	for ( Int_t ipmt=12; ipmt<24; ipmt++){
	183	if(time[ipmt] > 1) {
c41ceaac	184	if(time[ipmt]<besttimeA) {
	185	besttimeA=time[ipmt]; //timeA
	186	pmtBestA=ipmt;}
dc7ca31d	187	}
dc7ca31d	188	}
adf36b9d	189	if(besttimeA < 999999) {
8f620945	190	frecpoints->SetTimeBestA(Int_t(besttimeA *channelWidth - fdZonA/c));
adf36b9d	191	tr[1]=true;
	192	}
	193	if( besttimeC < 999999 ) {
8f620945	194	frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth - fdZonA/c));
adf36b9d	195	tr[2]=true;
adf36b9d	196	}
d0bcd1fb	197	AliDebug(10,Form(" besttimeA %f ch, besttimeC %f ch",besttimeA, besttimeC));
adf36b9d	198	if(besttimeA <999999 && besttimeC < 999999 ){
9b83615d	199	// timeDiff = (besttimeC - besttimeA)*channelWidth;
9b83615d	200	timeDiff = (besttimeA - besttimeC)*channelWidth;
adf36b9d	201	meanTime = (besttimeA + besttimeC)/2;// * channelWidth);
8f620945	202	timeclock = meanTime *channelWidth -fdZonA/c ;
adf36b9d	203	vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
	204	tr[0]=true;
	205	}
	206	frecpoints->SetVertex(vertex);
	207	frecpoints->SetMeanTime(meanTime);
	208	frecpoints->SetT0clock(timeclock);
	209	frecpoints->SetT0Trig(tr);
	210
8f620945	211	AliDebug(10,Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));
adf36b9d	212
	213	//online mean
	214	frecpoints->SetOnlineMean(Int_t(onlineMean));
b5a9f753	215	AliDebug(10,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm online mean %i timeclock %f ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
adf36b9d	216
	217
	218
	219
b95e8d87	220
dc7ca31d	221	clustersTree->Fill();
bd375212	222
	223	delete timeCFD;
	224	delete timeLED;
	225	delete chargeQT0;
	226	delete chargeQT1;
dc7ca31d	227	}
	228
	229
c41ceaac	230	//_______________________________________________________________________
	231
	232	void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
	233	{
94c27e4f	234	// T0 raw ->
539b9cb9	235	//
	236	// reference amplitude and time ref. point from reco param
	237
38cbfa7c	238	Float_t refAmp = GetRecoParam()->GetRefAmp();
	239	Int_t refPoint = 0;
	240	//Bad channel
	241	Int_t badpmt = GetRecoParam()->GetRefPoint();
	242
	243	Int_t low[110], high[110];
58bd3a16	244
e8ed1cd0	245	Int_t allData[110][5];
2e6a5ee0	246
e8ed1cd0	247	Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
adf36b9d	248	Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
	249	Float_t c = 29.9792458; // cm/ns
	250	Double32_t vertex = 9999999;
776de217	251	Int_t onlineMean=0;
9480f05f	252	// Float_t meanVertex = fParam->GetMeanVertex();
9480f05f	253	Float_t meanVertex = 0;
bce12dc5	254	for (Int_t i0=0; i0<105; i0++)
bce12dc5	255	{
38cbfa7c	256	for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
	257	low[i0] = Int_t (GetRecoParam()->GetLow(i0));
	258	high[i0] = Int_t (GetRecoParam()->GetHigh(i0));
bce12dc5	259	}
2e6a5ee0	260
adf36b9d	261	Double32_t besttimeA=9999999;
adf36b9d	262	Double32_t besttimeC=9999999;
bce12dc5	263	Int_t pmtBestA=99999;
bce12dc5	264	Int_t pmtBestC=99999;
29a60970	265
bce12dc5	266	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	267
	268	recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
2e6a5ee0	269
bce12dc5	270	AliDebug(10," before read data ");
bce12dc5	271	AliT0RawReader myrawreader(rawReader);
776de217	272
	273	UInt_t type =rawReader->GetType();
	274
bce12dc5	275	if (!myrawreader.Next())
	276	AliDebug(1,Form(" no raw data found!!"));
	277	else
	278	{
38cbfa7c	279	for (Int_t i=0; i<24; i++)
	280	{
	281	timeCFD[i]=0; timeLED[i]=0; chargeQT0[i]=0; chargeQT1[i]=0;
	282	}
	283	Int_t fBCID=Int_t (rawReader->GetBCID());
8f620945	284	Int_t trmbunch= myrawreader.GetTRMBunchID();
38cbfa7c	285	if(type == 7 ) { //only physics
8f620945	286	for (Int_t i=0; i<105; i++) {
bce12dc5	287	for (Int_t iHit=0; iHit<5; iHit++)
	288	{
	289	allData[i][iHit] = myrawreader.GetData(i,iHit);
	290	}
8f620945	291	}
	292	Int_t ref=0;
	293	if (refPoint>0)
	294	ref = allData[refPoint][0]-5000;
	295
	296	Float_t channelWidth = fParam->GetChannelWidth() ;
	297
	298	// Int_t meanT0 = fParam->GetMeanT0();
	299
38cbfa7c	300	for (Int_t in=0; in<12; in++)
	301	{
	302	for (Int_t iHit=0; iHit<5; iHit++)
	303	{
	304	if(allData[in+1][iHit] > low[in+1] &&
	305	allData[in+1][iHit] < high[in+1])
	306	{
	307	timeCFD[in] = allData[in+1][iHit] ;
	308	break;
	309	}
	310	}
	311	for (Int_t iHit=0; iHit<5; iHit++)
	312	{
	313	if(allData[in+12+1][iHit] > low[in+12+1] &&
	314	allData[in+1+12][iHit] < high[in+12+1])
	315	{
	316	timeLED[in] = allData[in+12+1][iHit] ;
	317	break;
	318	}
	319	}
	320	for (Int_t iHit=0; iHit<5; iHit++)
	321	{
	322	if(allData[in+1+56][iHit] > low[in+1+56] &&
	323	allData[in+1+56][iHit] < high[in+1+56])
	324	{
	325	timeCFD[in+12] = allData[in+56+1][iHit] ;
	326	break;
	327	}
	328	}
	329
	330	for (Int_t iHit=0; iHit<5; iHit++)
	331	{
	332	if(allData[in+1+68][iHit] > low[in+1+68] &&
	333	allData[in+1+68][iHit] < high[in+1+68])
	334	{
	335	timeLED[in+12] = allData[in+68+1][iHit] ;
	336	break;
	337	}
	338	}
	339	AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
	340	in, timeCFD[in],timeCFD[in+12],timeLED[in],
	341	timeLED[in+12]));
	342
	343	}
	344
	345
	346	for (Int_t iHit=0; iHit<5; iHit++)
8f620945	347	{
38cbfa7c	348	for (Int_t in=0; in<12; in++)
	349	{
	350	if(allData[2in+25][iHit] > low[2in+25] &&
	351	allData[2in+25][iHit] < high[2in+25])
	352	chargeQT0[in]=allData[2*in+25][iHit];
	353	if(allData[2in+26][iHit] > low[2in+26] &&
	354	allData[2in+26][iHit] < high[2in+26])
	355	chargeQT1[in]=allData[2*in+26][iHit];
	356
	357	}
	358
	359	for (Int_t in=12; in<24; in++)
	360	{
	361	if(allData[2in+57][iHit] > low[2in+57] &&
	362	allData[2in+57][iHit] < high[2in+57])
	363	chargeQT0[in]=allData[2*in+57][iHit];
	364	if(allData[2in+58][iHit] > low[2in+58] &&
	365	allData[2in+58][iHit] < high[2in+58])
	366	chargeQT1[in]=allData[2*in+58][iHit];
	367
	368	}
8f620945	369	}
8f620945	370
8f620945	371
8f620945	372	for (Int_t in=0; in<24; in++)
	373	AliDebug(10, Form(" readed Raw %i %i %i %i %i",
	374	in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
38cbfa7c	375	for (Int_t iHit=0; iHit<5; iHit++)
	376	{
	377	if(allData[49][iHit] > low[49] &&
	378	allData[49][iHit] < high[49])
	379	onlineMean = allData[49][iHit];
	380	}
73df58ab	381	Double32_t time[24], adc[24], noncalibtime[24];
8f620945	382	for (Int_t ipmt=0; ipmt<24; ipmt++) {
38cbfa7c	383	if(timeCFD[ipmt]>0 ){
bce12dc5	384	//for simulated data
d0bcd1fb	385	//for physics data
73df58ab	386	if(( chargeQT0[ipmt] - chargeQT1[ipmt])>0) {
541b42c4	387	adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
73df58ab	388	}
d0bcd1fb	389	else
d0bcd1fb	390	adc[ipmt] = 0;
38cbfa7c	391	time[ipmt] = fCalib-> WalkCorrection(Int_t (refAmp), ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
bce12dc5	392
d0bcd1fb	393	Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
8f620945	394	// time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
73df58ab	395	AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb	396	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
	397	Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	398	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
345f03db	399	AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %f\n ",ampMip,qtMip, adc[ipmt]));
73df58ab	400	//bad peak removing
8f620945	401	if(sl<550) {
73df58ab	402	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
73df58ab	403	// frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
345f03db	404	frecpoints->SetAmp(ipmt, Double32_t( qtMip)); //for cosmic &pp beam
345f03db	405	frecpoints->SetAmpLED(ipmt, Double32_t(ampMip));
73df58ab	406	noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
73df58ab	407	}
2e6a5ee0	408	}
	409	else {
	410	time[ipmt] = 0;
	411	adc[ipmt] = 0;
73df58ab	412	noncalibtime[ipmt] = 0;
2e6a5ee0	413	}
2e6a5ee0	414	}
73df58ab	415	fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;
2e6a5ee0	416	for (Int_t ipmt=0; ipmt<12; ipmt++){
38cbfa7c	417	if(time[ipmt] > 1 && ipmt != badpmt )
73df58ab	418	{
	419	if(time[ipmt]<besttimeC){
	420	besttimeC=time[ipmt]; //timeC
	421	pmtBestC=ipmt;
	422	}
2e6a5ee0	423	}
2e6a5ee0	424	}
73df58ab	425	for ( Int_t ipmt=12; ipmt<24; ipmt++)
73df58ab	426	{
38cbfa7c	427	if(time[ipmt] > 1 && ipmt != badpmt)
73df58ab	428	{
	429	if(time[ipmt]<besttimeA) {
	430	besttimeA=time[ipmt]; //timeA
	431	pmtBestA=ipmt;
	432	}
	433	}
2e6a5ee0	434	}
adf36b9d	435	if(besttimeA < 999999)
73df58ab	436	frecpoints->SetTimeBestA(besttimeA * channelWidth - 1000.fLatencyHPTDC + 1000.fLatencyL1A);
adf36b9d	437	if( besttimeC < 999999 )
73df58ab	438	frecpoints->SetTimeBestC(besttimeC * channelWidth - 1000.fLatencyHPTDC +1000.fLatencyL1C);
adf36b9d	439	AliDebug(10,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f ps",
d0bcd1fb	440	pmtBestA,besttimeA, pmtBestC, besttimeC));
adf36b9d	441	if(besttimeA <999999 && besttimeC < 999999 ){
73df58ab	442	timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
73df58ab	443	timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - 1000.fLatencyHPTDC + 1000.fLatencyL1;
adf36b9d	444	meanTime = (besttimeA+besttimeC-2.*Float_t(ref))/2.;
9480f05f	445	vertex = meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
adf36b9d	446	}
776de217	447	} //if phys event
adf36b9d	448	AliDebug(5,Form(" timeDiff %f #channel, meanTime %f #channel, TOFmean%f vertex %f cm meanVertex %f online mean %i \n",timeDiff, meanTime,timeclock, vertex,meanVertex, onlineMean));
776de217	449	frecpoints->SetT0clock(timeclock);
776de217	450	frecpoints->SetVertex(vertex);
adf36b9d	451	frecpoints->SetMeanTime(meanTime);
776de217	452	frecpoints->SetOnlineMean(Int_t(onlineMean));
adf36b9d	453	// Set triggers
	454
	455	Bool_t tr[5];
	456	Int_t trchan[5]= {50,51,52,55,56};
	457	for (Int_t i=0; i<5; i++) tr[i]=false;
	458	for (Int_t itr=0; itr<5; itr++) {
38cbfa7c	459	for (Int_t iHit=0; iHit<5; iHit++)
	460	{
	461	Int_t trr=trchan[itr];
	462	if(allData[trr][iHit] > low[trr] &&
	463	allData[trr][iHit] < high[trr]&&
	464	allData[trr][iHit] > 0) tr[itr]=true;
	465	}
adf36b9d	466	}
38cbfa7c	467	frecpoints->SetT0Trig(tr);
38cbfa7c	468	} // if (else )raw data
58bd3a16	469	recTree->Fill();
	470	if(frecpoints) delete frecpoints;
	471	}
adf36b9d	472
	473
	474	//____________________________________________________________
	475
	476	void AliT0Reconstructor::FillESD(TTree /digitsTree/, TTree clustersTree, AliESDEvent *pESD) const
	477	{
dc7ca31d	478
	479	/***************************************************
	480	Resonstruct digits to vertex position
	481	****************************************************/
	482
dc7ca31d	483	AliDebug(1,Form("Start FillESD T0"));
58bd3a16	484	Float_t channelWidth = fParam->GetChannelWidth() ;
f7c2c2fc	485	Float_t c = 0.0299792458; // cm/ps
adf36b9d	486	Float_t currentVertex=0, shift=0;
	487	Int_t ncont=0;
	488	const AliESDVertex* vertex = pESD->GetPrimaryVertex();
	489	if (!vertex) vertex = pESD->GetPrimaryVertexSPD();
	490	if (!vertex) vertex = pESD->GetPrimaryVertexTPC();
	491	if (!vertex) vertex = pESD->GetVertex();
	492
	493	if (vertex) {
	494	AliDebug(2, Form("Got %s (%s) from ESD: %f",
	495	vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
	496	currentVertex = vertex->GetZ();
	497
	498	ncont = vertex->GetNContributors();
	499	// cout<<" spdver "<<spdver<<" ncont "<<ncont<<endl;
58bd3a16	500	if(ncont>2 ) {
adf36b9d	501	shift = currentVertex/c;
	502	// cout<<" vertex shif "<<shift<<" vertex "<<spdver<<" IsFromVertexer3D "<<fverSPD->IsFromVertexer3D()<<endl;
	503	}
	504	}
d76c31f4	505	TTree *treeR = clustersTree;
dc7ca31d	506
73df58ab	507	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
73df58ab	508	if (!frecpoints) {
dc7ca31d	509	AliError("Reconstruct Fill ESD >> no recpoints found");
	510	return;
	511	}
	512
	513	AliDebug(1,Form("Start FillESD T0"));
	514	TBranch *brRec = treeR->GetBranch("T0");
	515	if (brRec) {
	516	brRec->SetAddress(&frecpoints);
	517	}else{
f16935f7	518	AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d	519	return;
dc7ca31d	520	}
73df58ab	521
	522	brRec->GetEntry(0);
	523	Double32_t amp[24], time[24], ampQTC[24], timecorr[24];
	524	Double32_t* tcorr;
	525	for(Int_t i=0; i<24; i++)
	526	amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
	527
	528	Double32_t timeClock[3];
	529	Double32_t zPosition = frecpoints -> GetVertex();
	530	Double32_t timeStart = frecpoints -> GetMeanTime();
	531	timeClock[0] = frecpoints -> GetT0clock() ;
	532	timeClock[1] = frecpoints -> GetBestTimeA() + shift;
	533	timeClock[2] = frecpoints -> GetBestTimeC() - shift;
	534	for ( Int_t i=0; i<24; i++) {
	535	time[i] = frecpoints -> GetTime(i); // ps to ns
	536	if ( time[i] >1) {
345f03db	537	ampQTC[i] = frecpoints -> GetAmp(i);
345f03db	538	amp[i] = frecpoints -> AmpLED(i);
dc7ca31d	539	}
73df58ab	540	}
	541	Int_t trig= frecpoints ->GetT0Trig();
	542	pESD->SetT0Trig(trig);
	543
	544	pESD->SetT0zVertex(zPosition); //vertex Z position
	545	pESD->SetT0(timeStart); // interaction time
	546	for(Int_t i=0; i<3; i++)
	547	pESD->SetT0TOF(i,timeClock[i]); // interaction time (ns)
	548	pESD->SetT0time(time); // best TOF on each PMT
345f03db	549	pESD->SetT0amplitude(ampQTC); // number of particles(MIPs) on each PMT
73df58ab	550
	551	AliDebug(1,Form("T0: Vertex %f (T0A+T0C)/2 %f #channels T0signal %f ns OrA %f ns OrC %f T0trig %i\n",zPosition, timeStart, timeClock[0], timeClock[1], timeClock[2], trig));
	552
	553	if (pESD) {
adf36b9d	554
73df58ab	555	AliESDfriend fr = (AliESDfriend)pESD->FindListObject("AliESDfriend");
	556	if (fr) {
	557	AliDebug(1, Form("Writing TZERO friend data to ESD tree"));
	558
	559	if (ncont>2) {
	560	tcorr = fESDTZEROfriend->GetT0timeCorr();
	561	for ( Int_t i=0; i<24; i++) {
	562	timecorr[i]=tcorr[i];
f7c2c2fc	563	if(i<12 && time[i]>1) timecorr[i] -= shift/channelWidth;
f7c2c2fc	564	if(i>11 && time[i]>1) timecorr[i] += shift/channelWidth;
58bd3a16	565	}
73df58ab	566	fESDTZEROfriend->SetT0timeCorr( timecorr) ;
	567	fESDTZEROfriend->SetT0ampLEDminCFD(amp);
	568	fESDTZEROfriend->SetT0ampQTC(ampQTC);
	569	fr->SetTZEROfriend(fESDTZEROfriend);
	570	}//
58bd3a16	571	}
58bd3a16	572	}
73df58ab	573
	574
	575
dc7ca31d	576	} // vertex in 3 sigma
	577
	578
	579
	580
	581
	582