[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 <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()
{
  //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) ; 	
  }

  
  // fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
  //fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));


  fCalib = new AliT0Calibrator();

}

//_____________________________________________________________________________
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);

  for (Int_t i=0; i<5; i++) {
    printf(" T0 trigers %i ",tr[i]);
  }
    printf(" \n ");

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