[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"));


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

}

//_____________________________________________________________________________
void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
  
{
  // T0 digits reconstruction
  Int_t refAmp = 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, 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 %i\n ",ampMip,qtMip, adc[ipmt]));
      
      frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
      frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam 
      frecpoints->SetAmpLED(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));
    tr[1]=true;
  }
  if( besttimeC < 999999 ) {
    frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth));
    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 ;
    vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
    tr[0]=true; 
  }
  frecpoints->SetVertex(vertex);
  frecpoints->SetMeanTime(meanTime);
  frecpoints->SetT0clock(timeclock);
  frecpoints->SetT0Trig(tr);

  AliInfo(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 %i #channel,  meanTime %i #channel, vertex %f cm online mean %i timeclock %i 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

  Int_t refAmp = GetRecoParam()->GetRefAmp();
  Int_t refPoint = GetRecoParam()->GetRefPoint();

  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; 	
    }
   
  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
    {  
      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++)  
	    {
	      timeCFD[in] = allData[in+1][0] ;
	      timeCFD[in+12] = allData[in+56+1][0] ;
	      timeLED[in] = allData[in+12+1][0] ;
	      timeLED[in+12] = allData[in+68+1][0] ;
	      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 in=0; in<12;  in++)
	    {
	      chargeQT0[in]=allData[2*in+25][0];
	      chargeQT1[in]=allData[2*in+26][0];
	    }
	  
	   for (Int_t in=12; in<24;  in++)
	     {
	       chargeQT0[in]=allData[2*in+57][0];
	       chargeQT1[in]=allData[2*in+58][0];
	     }
	   
	   //	 } //cosmic with physics event
       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]));
        onlineMean = allData[49][0];       
       
	Double32_t time[24], adc[24],  noncalibtime[24];
       for (Int_t ipmt=0; ipmt<24; ipmt++) {
	 if(timeCFD[ipmt]>0 && timeLED[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(refAmp, ipmt, 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(0,Form("  Amlitude in MIPS LED %f ; QTC %f;  in channels %i\n ",ampMip,qtMip, adc[ipmt]));
	   //bad peak removing
	   if(sl<540) {
	     frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	     // frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
	     frecpoints->SetAmpLED(ipmt, Double32_t( qtMip)); //for cosmic &pp beam 
	     frecpoints->SetAmp(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 
	    && (timeLED[ipmt] - timeCFD[ipmt])<550 )
	   {
	     if(time[ipmt]<besttimeC){
	       besttimeC=time[ipmt]; //timeC
	       pmtBestC=ipmt;
	     }
	   }
       }
       for ( Int_t ipmt=12; ipmt<24; ipmt++)
	 {
	   if(time[ipmt] > 1 
	      && (timeLED[ipmt] - timeCFD[ipmt])<550 )
	     {
	       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++) {
	if(allData[trchan[itr]][0]>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 = 29.9792458; // cm/ns
  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) {
      amp[i] = frecpoints -> GetAmp(i);
      ampQTC[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(amp);     // 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"));
	82	//fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
	83
539b9cb9	84
12e9daf9	85	fCalib = new AliT0Calibrator();
73df58ab	86	fESDTZEROfriend = new AliESDTZEROfriend();
12e9daf9	87
dc7ca31d	88	}
c41ceaac	89
c41ceaac	90	//_____________________________________________________________________________
dc7ca31d	91	void AliT0Reconstructor::Reconstruct(TTreedigitsTree, TTreeclustersTree) const
94c27e4f	92
dc7ca31d	93	{
94c27e4f	94	// T0 digits reconstruction
c555f418	95	Int_t refAmp = GetRecoParam()->GetRefAmp();
776de217	96
c41ceaac	97	TArrayI * timeCFD = new TArrayI(24);
	98	TArrayI * timeLED = new TArrayI(24);
	99	TArrayI * chargeQT0 = new TArrayI(24);
	100	TArrayI * chargeQT1 = new TArrayI(24);
74adb36a	101
d0bcd1fb	102
8955c6b4	103	Float_t channelWidth = fParam->GetChannelWidth() ;
b95e8d87	104	Float_t meanVertex = fParam->GetMeanVertex();
776de217	105	Float_t c = 0.0299792; // cm/ps
adf36b9d	106	Double32_t vertex = 9999999;
adf36b9d	107	Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
776de217	108
94c27e4f	109
dc7ca31d	110	AliDebug(1,Form("Start DIGITS reconstruction "));
94c27e4f	111
2e6a5ee0	112
d0bcd1fb	113	TBranch *brDigits=digitsTree->GetBranch("T0");
e0bba6cc	114	AliT0digit *fDigits = new AliT0digit() ;
dc7ca31d	115	if (brDigits) {
	116	brDigits->SetAddress(&fDigits);
	117	}else{
f16935f7	118	AliError(Form("EXEC Branch T0 digits not found"));
f16935f7	119	return;
dc7ca31d	120	}
e0bba6cc	121
c41ceaac	122	digitsTree->GetEvent(0);
	123	digitsTree->GetEntry(0);
	124	brDigits->GetEntry(0);
	125	fDigits->GetTimeCFD(*timeCFD);
	126	fDigits->GetTimeLED(*timeLED);
	127	fDigits->GetQT0(*chargeQT0);
	128	fDigits->GetQT1(*chargeQT1);
446d6ec4	129	Int_t onlineMean = fDigits->MeanTime();
c883fdf2	130
adf36b9d	131	Bool_t tr[5];
adf36b9d	132	for (Int_t i=0; i<5; i++) tr[i]=false;
c41ceaac	133
adf36b9d	134	Double32_t besttimeA=999999;
adf36b9d	135	Double32_t besttimeC=999999;
c41ceaac	136	Int_t pmtBestA=99999;
c41ceaac	137	Int_t pmtBestC=99999;
dc7ca31d	138
94c27e4f	139	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	140	clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
	141
b95e8d87	142	Float_t time[24], adc[24];
dc7ca31d	143	for (Int_t ipmt=0; ipmt<24; ipmt++) {
c41ceaac	144	if(timeCFD->At(ipmt)>0 ){
d0bcd1fb	145	if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)
d0bcd1fb	146	adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
a7027400	147	else
d0bcd1fb	148	adc[ipmt] = 0;
d0bcd1fb	149
541b42c4	150	// time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD->At(ipmt)) ;
d0bcd1fb	151
d0bcd1fb	152	Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
541b42c4	153	time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD->At(ipmt) ) ;
d0bcd1fb	154	AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb	155	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
776de217	156
d0bcd1fb	157	Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	158	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
	159	AliDebug(10,Form(" Amlitude in MIPS LED %f , QTC %f in channels %i\n ",ampMip,qtMip, adc[ipmt]));
	160
	161	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	162	frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam
	163	frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
	164
dc7ca31d	165	}
	166	else {
	167	time[ipmt] = 0;
	168	adc[ipmt] = 0;
	169	}
	170	}
94c27e4f	171
dc7ca31d	172	for (Int_t ipmt=0; ipmt<12; ipmt++){
dc7ca31d	173	if(time[ipmt] > 1 ) {
c41ceaac	174	if(time[ipmt]<besttimeC){
	175	besttimeC=time[ipmt]; //timeC
	176	pmtBestC=ipmt;
dc7ca31d	177	}
	178	}
	179	}
	180	for ( Int_t ipmt=12; ipmt<24; ipmt++){
	181	if(time[ipmt] > 1) {
c41ceaac	182	if(time[ipmt]<besttimeA) {
	183	besttimeA=time[ipmt]; //timeA
	184	pmtBestA=ipmt;}
dc7ca31d	185	}
dc7ca31d	186	}
adf36b9d	187	if(besttimeA < 999999) {
	188	frecpoints->SetTimeBestA(Int_t(besttimeA *channelWidth));
	189	tr[1]=true;
	190	}
	191	if( besttimeC < 999999 ) {
	192	frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth));
	193	tr[2]=true;
	194	}
d0bcd1fb	195	AliDebug(10,Form(" besttimeA %f ch, besttimeC %f ch",besttimeA, besttimeC));
adf36b9d	196	if(besttimeA <999999 && besttimeC < 999999 ){
9b83615d	197	// timeDiff = (besttimeC - besttimeA)*channelWidth;
9b83615d	198	timeDiff = (besttimeA - besttimeC)*channelWidth;
adf36b9d	199	meanTime = (besttimeA + besttimeC)/2;// * channelWidth);
	200	timeclock = meanTime *channelWidth ;
	201	vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
	202	tr[0]=true;
	203	}
	204	frecpoints->SetVertex(vertex);
	205	frecpoints->SetMeanTime(meanTime);
	206	frecpoints->SetT0clock(timeclock);
	207	frecpoints->SetT0Trig(tr);
	208
8f8d0732	209	AliInfo(Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));
adf36b9d	210
	211	//online mean
	212	frecpoints->SetOnlineMean(Int_t(onlineMean));
	213	AliDebug(10,Form(" timeDiff %i #channel, meanTime %i #channel, vertex %f cm online mean %i timeclock %i ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
	214
	215
	216
	217
b95e8d87	218
dc7ca31d	219	clustersTree->Fill();
bd375212	220
	221	delete timeCFD;
	222	delete timeLED;
	223	delete chargeQT0;
	224	delete chargeQT1;
dc7ca31d	225	}
	226
	227
c41ceaac	228	//_______________________________________________________________________
	229
	230	void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
	231	{
94c27e4f	232	// T0 raw ->
539b9cb9	233	//
	234	// reference amplitude and time ref. point from reco param
	235
c555f418	236	Int_t refAmp = GetRecoParam()->GetRefAmp();
c555f418	237	Int_t refPoint = GetRecoParam()->GetRefPoint();
58bd3a16	238
e8ed1cd0	239	Int_t allData[110][5];
2e6a5ee0	240
e8ed1cd0	241	Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
adf36b9d	242	Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
	243	Float_t c = 29.9792458; // cm/ns
	244	Double32_t vertex = 9999999;
776de217	245	Int_t onlineMean=0;
9480f05f	246	// Float_t meanVertex = fParam->GetMeanVertex();
9480f05f	247	Float_t meanVertex = 0;
bce12dc5	248	for (Int_t i0=0; i0<105; i0++)
	249	{
	250	for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
	251	}
2e6a5ee0	252
adf36b9d	253	Double32_t besttimeA=9999999;
adf36b9d	254	Double32_t besttimeC=9999999;
bce12dc5	255	Int_t pmtBestA=99999;
bce12dc5	256	Int_t pmtBestC=99999;
29a60970	257
bce12dc5	258	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	259
	260	recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
2e6a5ee0	261
bce12dc5	262	AliDebug(10," before read data ");
bce12dc5	263	AliT0RawReader myrawreader(rawReader);
776de217	264
	265	UInt_t type =rawReader->GetType();
	266
bce12dc5	267	if (!myrawreader.Next())
	268	AliDebug(1,Form(" no raw data found!!"));
	269	else
	270	{
776de217	271	if(type == 7) { //only physics
bce12dc5	272	for (Int_t i=0; i<105; i++) {
	273	for (Int_t iHit=0; iHit<5; iHit++)
	274	{
	275	allData[i][iHit] = myrawreader.GetData(i,iHit);
	276	}
	277	}
541b42c4	278	Int_t ref=0;
541b42c4	279	if (refPoint>0)
adf36b9d	280	ref = allData[refPoint][0]-5000;
541b42c4	281
bce12dc5	282	Float_t channelWidth = fParam->GetChannelWidth() ;
	283
	284	// Int_t meanT0 = fParam->GetMeanT0();
d0bcd1fb	285
bce12dc5	286
	287	for (Int_t in=0; in<12; in++)
	288	{
	289	timeCFD[in] = allData[in+1][0] ;
	290	timeCFD[in+12] = allData[in+56+1][0] ;
	291	timeLED[in] = allData[in+12+1][0] ;
	292	timeLED[in+12] = allData[in+68+1][0] ;
d0bcd1fb	293	AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
	294	in, timeCFD[in],timeCFD[in+12],timeLED[in],
	295	timeLED[in+12]));
bce12dc5	296	}
	297
	298	for (Int_t in=0; in<12; in++)
	299	{
d0bcd1fb	300	chargeQT0[in]=allData[2*in+25][0];
d0bcd1fb	301	chargeQT1[in]=allData[2*in+26][0];
bce12dc5	302	}
	303
	304	for (Int_t in=12; in<24; in++)
	305	{
d0bcd1fb	306	chargeQT0[in]=allData[2*in+57][0];
d0bcd1fb	307	chargeQT1[in]=allData[2*in+58][0];
bce12dc5	308	}
bce12dc5	309
d0bcd1fb	310	// } //cosmic with physics event
2e6a5ee0	311	for (Int_t in=0; in<24; in++)
c883fdf2	312	AliDebug(10, Form(" readed Raw %i %i %i %i %i",
c883fdf2	313	in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
776de217	314	onlineMean = allData[49][0];
2e6a5ee0	315
73df58ab	316	Double32_t time[24], adc[24], noncalibtime[24];
2e6a5ee0	317	for (Int_t ipmt=0; ipmt<24; ipmt++) {
2e6a5ee0	318	if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
bce12dc5	319	//for simulated data
d0bcd1fb	320	//for physics data
73df58ab	321	if(( chargeQT0[ipmt] - chargeQT1[ipmt])>0) {
541b42c4	322	adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
73df58ab	323	}
d0bcd1fb	324	else
d0bcd1fb	325	adc[ipmt] = 0;
bce12dc5	326
539b9cb9	327
40043664	328	// time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
bce12dc5	329
d0bcd1fb	330	Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
40043664	331	time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
73df58ab	332	AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb	333	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
	334	Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	335	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
73df58ab	336	AliDebug(0,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %i\n ",ampMip,qtMip, adc[ipmt]));
	337	//bad peak removing
	338	if(sl<540) {
	339	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	340	// frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
	341	frecpoints->SetAmpLED(ipmt, Double32_t( qtMip)); //for cosmic &pp beam
	342	frecpoints->SetAmp(ipmt, Double32_t(ampMip));
	343	noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
	344	}
2e6a5ee0	345	}
	346	else {
	347	time[ipmt] = 0;
	348	adc[ipmt] = 0;
73df58ab	349	noncalibtime[ipmt] = 0;
2e6a5ee0	350	}
2e6a5ee0	351	}
73df58ab	352	fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;
2e6a5ee0	353	for (Int_t ipmt=0; ipmt<12; ipmt++){
73df58ab	354	if(time[ipmt] > 1
	355	&& (timeLED[ipmt] - timeCFD[ipmt])<550 )
	356	{
	357	if(time[ipmt]<besttimeC){
	358	besttimeC=time[ipmt]; //timeC
	359	pmtBestC=ipmt;
	360	}
2e6a5ee0	361	}
2e6a5ee0	362	}
73df58ab	363	for ( Int_t ipmt=12; ipmt<24; ipmt++)
	364	{
	365	if(time[ipmt] > 1
	366	&& (timeLED[ipmt] - timeCFD[ipmt])<550 )
	367	{
	368	if(time[ipmt]<besttimeA) {
	369	besttimeA=time[ipmt]; //timeA
	370	pmtBestA=ipmt;
	371	}
	372	}
2e6a5ee0	373	}
adf36b9d	374	if(besttimeA < 999999)
73df58ab	375	frecpoints->SetTimeBestA(besttimeA * channelWidth - 1000.fLatencyHPTDC + 1000.fLatencyL1A);
adf36b9d	376	if( besttimeC < 999999 )
73df58ab	377	frecpoints->SetTimeBestC(besttimeC * channelWidth - 1000.fLatencyHPTDC +1000.fLatencyL1C);
adf36b9d	378	AliDebug(10,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f ps",
d0bcd1fb	379	pmtBestA,besttimeA, pmtBestC, besttimeC));
adf36b9d	380	if(besttimeA <999999 && besttimeC < 999999 ){
73df58ab	381	timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
73df58ab	382	timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - 1000.fLatencyHPTDC + 1000.fLatencyL1;
adf36b9d	383	meanTime = (besttimeA+besttimeC-2.*Float_t(ref))/2.;
9480f05f	384	vertex = meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
adf36b9d	385	}
776de217	386	} //if phys event
adf36b9d	387	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	388	frecpoints->SetT0clock(timeclock);
776de217	389	frecpoints->SetVertex(vertex);
adf36b9d	390	frecpoints->SetMeanTime(meanTime);
776de217	391	frecpoints->SetOnlineMean(Int_t(onlineMean));
adf36b9d	392	// Set triggers
	393
	394	Bool_t tr[5];
	395	Int_t trchan[5]= {50,51,52,55,56};
	396	for (Int_t i=0; i<5; i++) tr[i]=false;
	397	for (Int_t itr=0; itr<5; itr++) {
	398	if(allData[trchan[itr]][0]>0) tr[itr]=true;
	399	frecpoints->SetT0Trig(tr);
	400	}
58bd3a16	401	} // if (else )raw data
	402	recTree->Fill();
	403	if(frecpoints) delete frecpoints;
	404	}
adf36b9d	405
	406
	407	//____________________________________________________________
	408
	409	void AliT0Reconstructor::FillESD(TTree /digitsTree/, TTree clustersTree, AliESDEvent *pESD) const
	410	{
dc7ca31d	411
	412	/***************************************************
	413	Resonstruct digits to vertex position
	414	****************************************************/
	415
dc7ca31d	416	AliDebug(1,Form("Start FillESD T0"));
58bd3a16	417	Float_t channelWidth = fParam->GetChannelWidth() ;
adf36b9d	418	Float_t c = 29.9792458; // cm/ns
	419	Float_t currentVertex=0, shift=0;
	420	Int_t ncont=0;
	421	const AliESDVertex* vertex = pESD->GetPrimaryVertex();
	422	if (!vertex) vertex = pESD->GetPrimaryVertexSPD();
	423	if (!vertex) vertex = pESD->GetPrimaryVertexTPC();
	424	if (!vertex) vertex = pESD->GetVertex();
	425
	426	if (vertex) {
	427	AliDebug(2, Form("Got %s (%s) from ESD: %f",
	428	vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
	429	currentVertex = vertex->GetZ();
	430
	431	ncont = vertex->GetNContributors();
	432	// cout<<" spdver "<<spdver<<" ncont "<<ncont<<endl;
58bd3a16	433	if(ncont>2 ) {
adf36b9d	434	shift = currentVertex/c;
	435	// cout<<" vertex shif "<<shift<<" vertex "<<spdver<<" IsFromVertexer3D "<<fverSPD->IsFromVertexer3D()<<endl;
	436	}
	437	}
d76c31f4	438	TTree *treeR = clustersTree;
dc7ca31d	439
73df58ab	440	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
73df58ab	441	if (!frecpoints) {
dc7ca31d	442	AliError("Reconstruct Fill ESD >> no recpoints found");
	443	return;
	444	}
	445
	446	AliDebug(1,Form("Start FillESD T0"));
	447	TBranch *brRec = treeR->GetBranch("T0");
	448	if (brRec) {
	449	brRec->SetAddress(&frecpoints);
	450	}else{
f16935f7	451	AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d	452	return;
dc7ca31d	453	}
73df58ab	454
	455	brRec->GetEntry(0);
	456	Double32_t amp[24], time[24], ampQTC[24], timecorr[24];
	457	Double32_t* tcorr;
	458	for(Int_t i=0; i<24; i++)
	459	amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
	460
	461	Double32_t timeClock[3];
	462	Double32_t zPosition = frecpoints -> GetVertex();
	463	Double32_t timeStart = frecpoints -> GetMeanTime();
	464	timeClock[0] = frecpoints -> GetT0clock() ;
	465	timeClock[1] = frecpoints -> GetBestTimeA() + shift;
	466	timeClock[2] = frecpoints -> GetBestTimeC() - shift;
	467	for ( Int_t i=0; i<24; i++) {
	468	time[i] = frecpoints -> GetTime(i); // ps to ns
	469	if ( time[i] >1) {
dc7ca31d	470	amp[i] = frecpoints -> GetAmp(i);
58bd3a16	471	ampQTC[i] = frecpoints -> AmpLED(i);
dc7ca31d	472	}
73df58ab	473	}
	474	Int_t trig= frecpoints ->GetT0Trig();
	475	pESD->SetT0Trig(trig);
	476
	477	pESD->SetT0zVertex(zPosition); //vertex Z position
	478	pESD->SetT0(timeStart); // interaction time
	479	for(Int_t i=0; i<3; i++)
	480	pESD->SetT0TOF(i,timeClock[i]); // interaction time (ns)
	481	pESD->SetT0time(time); // best TOF on each PMT
	482	pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT
	483
	484	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));
	485
	486	if (pESD) {
adf36b9d	487
73df58ab	488	AliESDfriend fr = (AliESDfriend)pESD->FindListObject("AliESDfriend");
	489	if (fr) {
	490	AliDebug(1, Form("Writing TZERO friend data to ESD tree"));
	491
	492	if (ncont>2) {
	493	tcorr = fESDTZEROfriend->GetT0timeCorr();
	494	for ( Int_t i=0; i<24; i++) {
	495	timecorr[i]=tcorr[i];
	496	if(i<12 && time[i]>1) timecorr[i] -= shift*channelWidth;
	497	if(i>11 && time[i]>1) timecorr[i] += shift*channelWidth;
58bd3a16	498	}
73df58ab	499	fESDTZEROfriend->SetT0timeCorr( timecorr) ;
	500	fESDTZEROfriend->SetT0ampLEDminCFD(amp);
	501	fESDTZEROfriend->SetT0ampQTC(ampQTC);
	502	fr->SetTZEROfriend(fESDTZEROfriend);
	503	}//
58bd3a16	504	}
58bd3a16	505	}
73df58ab	506
	507
	508
dc7ca31d	509	} // vertex in 3 sigma
	510
	511
	512
	513
	514
	515