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

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

}

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

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