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

}
//____________________________________________________________________

AliT0Reconstructor::AliT0Reconstructor(const AliT0Reconstructor &r):
  AliReconstructor(r),
					     fdZonA(0),
	                                     fdZonC(0),
					     fZposition(0),
					     fParam(NULL),
                                             fAmpLEDrec(0),
					     fQTC(0),
					     fAmpLED(0),
  					     fCalib()

 {
  //
  // AliT0Reconstructor copy constructor
  //

  ((AliT0Reconstructor &) r).Copy(*this);

}

//_____________________________________________________________________________
AliT0Reconstructor &AliT0Reconstructor::operator=(const AliT0Reconstructor &r)
{
  //
  // Assignment operator
  //
  if (this != &r) ((AliT0Reconstructor &) r).Copy(*this);
  return *this;

}

//_____________________________________________________________________________

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
  Float_t vertex = 9999999;
  Int_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();
  Int_t ref =  fDigits->RefPoint();

  
  Float_t besttimeA=999999;
  Float_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( ipmt, Int_t(sl), timeCFD[ipmt],"cosmic" ) ;
      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));
  if( besttimeC != 999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
  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 = Float_t((besttimeA + besttimeC)/2);// * channelWidth); 
    timeclock = Float_t(meanTime - ref);
    vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2; 
   }
    frecpoints->SetVertex(vertex);
    frecpoints->SetMeanTime(meanTime);
    frecpoints->SetT0clock(timeclock);
    //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];
  Int_t timeDiff=999999, meanTime=999999, timeclock=999999;
  Float_t c = 0.0299792458; // cm/ps
  Float_t vertex = 9999999;
  Int_t onlineMean=0;
  for (Int_t i0=0; i0<105; i0++)
    {
      for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0; 	
    }
   
  Float_t besttimeA=9999999;
  Float_t besttimeC=9999999;
  Int_t pmtBestA=99999;
  Int_t pmtBestC=99999;
  Float_t meanVertex = fParam->GetMeanVertex();
   
  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 = 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];       
       
       Float_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] = chargeQT1[ipmt] - chargeQT0[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( 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->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 !=9999999)  frecpoints->SetTimeBestA(Int_t(besttimeA));
       if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
       AliDebug(5,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f #channel",
		       pmtBestA,besttimeA, pmtBestC,  besttimeC));
        if(besttimeA <9999999 && besttimeC < 9999999 ){
	 timeDiff = ( besttimeA - besttimeC) *channelWidth;
	 meanTime =  Float_t((besttimeA + besttimeC)/2.);  
	 timeclock = Float_t(meanTime - ref);
	 vertex =  meanVertex - c*(timeDiff)/2.; //+ (fdZonA - fdZonC)/2; 
       }
      }  //if phys event       
      frecpoints->SetT0clock(timeclock);
      frecpoints->SetVertex(vertex);
      frecpoints->SetMeanTime(Int_t(meanTime));
      frecpoints->SetOnlineMean(Int_t(onlineMean));
      AliDebug(5,Form("  timeDiff %f #channel,  meanTime %f #channel, vertex %f cm meanVertex %f online mean %i ",timeDiff, meanTime,vertex,meanVertex, onlineMean));
	 
       
    } // 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"));

  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);
    Float_t amp[24], time[24];
    Float_t  zPosition = frecpoints -> GetVertex();
    Float_t timeStart = frecpoints -> GetMeanTime() ;
    Float_t timeClock = frecpoints -> GetT0clock() ;
     for ( Int_t i=0; i<24; i++) {
      time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
      amp[i] = frecpoints -> GetAmp(i);
    }
    pESD->SetT0zVertex(zPosition); //vertex Z position 
    pESD->SetT0(timeStart);        // interaction time 
    pESD->SetT0clock(timeClock);   // interaction time with ref.point(spectrum) 
    pESD->SetT0time(time);         // best TOF on each PMT 
    pESD->SetT0amplitude(amp);     // number of particles(MIPs) on each PMT
 
    AliDebug(1,Form(" Z position %f cm,  T0  %f ps",zPosition , timeStart));

} // 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
f16935f7	66	fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
f16935f7	67	fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
539b9cb9	68
12e9daf9	69	fCalib = new AliT0Calibrator();
12e9daf9	70
c41ceaac	71	}
c41ceaac	72	//____________________________________________________________________
e0bba6cc	73
c41ceaac	74	AliT0Reconstructor::AliT0Reconstructor(const AliT0Reconstructor &r):
d76c31f4	75	AliReconstructor(r),
f16935f7	76	fdZonA(0),
539b9cb9	77	fdZonC(0),
f16935f7	78	fZposition(0),
f16935f7	79	fParam(NULL),
539b9cb9	80	fAmpLEDrec(0),
c883fdf2	81	fQTC(0),
c883fdf2	82	fAmpLED(0),
539b9cb9	83	fCalib()
2e6a5ee0	84
f16935f7	85	{
c41ceaac	86	//
	87	// AliT0Reconstructor copy constructor
	88	//
e0bba6cc	89
c41ceaac	90	((AliT0Reconstructor &) r).Copy(*this);
e0bba6cc	91
	92	}
	93
c41ceaac	94	//_____________________________________________________________________________
c41ceaac	95	AliT0Reconstructor &AliT0Reconstructor::operator=(const AliT0Reconstructor &r)
dc7ca31d	96	{
c41ceaac	97	//
	98	// Assignment operator
	99	//
c41ceaac	100	if (this != &r) ((AliT0Reconstructor &) r).Copy(*this);
c41ceaac	101	return *this;
e0bba6cc	102
dc7ca31d	103	}
c41ceaac	104
	105	//_____________________________________________________________________________
	106
dc7ca31d	107	void AliT0Reconstructor::Reconstruct(TTreedigitsTree, TTreeclustersTree) const
94c27e4f	108
dc7ca31d	109	{
94c27e4f	110	// T0 digits reconstruction
c555f418	111	Int_t refAmp = GetRecoParam()->GetRefAmp();
776de217	112
c41ceaac	113	TArrayI * timeCFD = new TArrayI(24);
	114	TArrayI * timeLED = new TArrayI(24);
	115	TArrayI * chargeQT0 = new TArrayI(24);
	116	TArrayI * chargeQT1 = new TArrayI(24);
74adb36a	117
d0bcd1fb	118
8955c6b4	119	Float_t channelWidth = fParam->GetChannelWidth() ;
b95e8d87	120	Float_t meanVertex = fParam->GetMeanVertex();
776de217	121	Float_t c = 0.0299792; // cm/ps
	122	Float_t vertex = 9999999;
	123	Int_t timeDiff=999999, meanTime=999999, timeclock=999999;
	124
94c27e4f	125
dc7ca31d	126	AliDebug(1,Form("Start DIGITS reconstruction "));
94c27e4f	127
2e6a5ee0	128
d0bcd1fb	129	TBranch *brDigits=digitsTree->GetBranch("T0");
e0bba6cc	130	AliT0digit *fDigits = new AliT0digit() ;
dc7ca31d	131	if (brDigits) {
	132	brDigits->SetAddress(&fDigits);
	133	}else{
f16935f7	134	AliError(Form("EXEC Branch T0 digits not found"));
f16935f7	135	return;
dc7ca31d	136	}
e0bba6cc	137
c41ceaac	138	digitsTree->GetEvent(0);
	139	digitsTree->GetEntry(0);
	140	brDigits->GetEntry(0);
	141	fDigits->GetTimeCFD(*timeCFD);
	142	fDigits->GetTimeLED(*timeLED);
	143	fDigits->GetQT0(*chargeQT0);
	144	fDigits->GetQT1(*chargeQT1);
446d6ec4	145	Int_t onlineMean = fDigits->MeanTime();
776de217	146	Int_t ref = fDigits->RefPoint();
c883fdf2	147
c41ceaac	148
	149	Float_t besttimeA=999999;
	150	Float_t besttimeC=999999;
	151	Int_t pmtBestA=99999;
	152	Int_t pmtBestC=99999;
dc7ca31d	153
94c27e4f	154	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	155	clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
	156
b95e8d87	157	Float_t time[24], adc[24];
dc7ca31d	158	for (Int_t ipmt=0; ipmt<24; ipmt++) {
c41ceaac	159	if(timeCFD->At(ipmt)>0 ){
d0bcd1fb	160	if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)
d0bcd1fb	161	adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
a7027400	162	else
d0bcd1fb	163	adc[ipmt] = 0;
d0bcd1fb	164
c555f418	165	time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD->At(ipmt)) ;
d0bcd1fb	166
	167	Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
	168	// time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt],"cosmic" ) ;
	169	AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
	170	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
776de217	171
d0bcd1fb	172	Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	173	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
	174	AliDebug(10,Form(" Amlitude in MIPS LED %f , QTC %f in channels %i\n ",ampMip,qtMip, adc[ipmt]));
	175
	176	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	177	frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam
	178	frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
	179
dc7ca31d	180	}
	181	else {
	182	time[ipmt] = 0;
	183	adc[ipmt] = 0;
	184	}
	185	}
94c27e4f	186
dc7ca31d	187	for (Int_t ipmt=0; ipmt<12; ipmt++){
dc7ca31d	188	if(time[ipmt] > 1 ) {
c41ceaac	189	if(time[ipmt]<besttimeC){
	190	besttimeC=time[ipmt]; //timeC
	191	pmtBestC=ipmt;
dc7ca31d	192	}
	193	}
	194	}
	195	for ( Int_t ipmt=12; ipmt<24; ipmt++){
	196	if(time[ipmt] > 1) {
c41ceaac	197	if(time[ipmt]<besttimeA) {
	198	besttimeA=time[ipmt]; //timeA
	199	pmtBestA=ipmt;}
dc7ca31d	200	}
dc7ca31d	201	}
c41ceaac	202	if(besttimeA !=999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
c41ceaac	203	if( besttimeC != 999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
d0bcd1fb	204	AliDebug(10,Form(" besttimeA %f ch, besttimeC %f ch",besttimeA, besttimeC));
776de217	205	if(besttimeA !=999999 && besttimeC != 999999 ){
9b83615d	206	// timeDiff = (besttimeC - besttimeA)*channelWidth;
9b83615d	207	timeDiff = (besttimeA - besttimeC)*channelWidth;
d0bcd1fb	208	meanTime = Float_t((besttimeA + besttimeC)/2);// * channelWidth);
776de217	209	timeclock = Float_t(meanTime - ref);
d0bcd1fb	210	vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
776de217	211	}
dc7ca31d	212	frecpoints->SetVertex(vertex);
776de217	213	frecpoints->SetMeanTime(meanTime);
776de217	214	frecpoints->SetT0clock(timeclock);
446d6ec4	215	//online mean
d0bcd1fb	216	frecpoints->SetOnlineMean(Int_t(onlineMean));
776de217	217	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));
dc7ca31d	218
776de217	219	// }
b95e8d87	220
dc7ca31d	221	clustersTree->Fill();
bd375212	222
	223	delete timeCFD;
	224	delete timeLED;
	225	delete chargeQT0;
	226	delete chargeQT1;
dc7ca31d	227	}
	228
	229
c41ceaac	230	//_______________________________________________________________________
	231
	232	void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
	233	{
94c27e4f	234	// T0 raw ->
539b9cb9	235	//
	236	// reference amplitude and time ref. point from reco param
	237
c555f418	238	Int_t refAmp = GetRecoParam()->GetRefAmp();
c555f418	239	Int_t refPoint = GetRecoParam()->GetRefPoint();
539b9cb9	240
e8ed1cd0	241	Int_t allData[110][5];
2e6a5ee0	242
e8ed1cd0	243	Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
776de217	244	Int_t timeDiff=999999, meanTime=999999, timeclock=999999;
	245	Float_t c = 0.0299792458; // cm/ps
	246	Float_t vertex = 9999999;
	247	Int_t onlineMean=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
bce12dc5	253	Float_t besttimeA=9999999;
	254	Float_t besttimeC=9999999;
	255	Int_t pmtBestA=99999;
	256	Int_t pmtBestC=99999;
b95e8d87	257	Float_t meanVertex = fParam->GetMeanVertex();
29a60970	258
bce12dc5	259	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	260
	261	recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
2e6a5ee0	262
bce12dc5	263	AliDebug(10," before read data ");
bce12dc5	264	AliT0RawReader myrawreader(rawReader);
776de217	265
	266	UInt_t type =rawReader->GetType();
	267
bce12dc5	268	if (!myrawreader.Next())
	269	AliDebug(1,Form(" no raw data found!!"));
	270	else
	271	{
776de217	272	if(type == 7) { //only physics
bce12dc5	273	for (Int_t i=0; i<105; i++) {
	274	for (Int_t iHit=0; iHit<5; iHit++)
	275	{
	276	allData[i][iHit] = myrawreader.GetData(i,iHit);
	277	}
	278	}
	279
c555f418	280	Int_t ref = allData[refPoint][0]-5000.;
bce12dc5	281	Float_t channelWidth = fParam->GetChannelWidth() ;
	282
	283	// Int_t meanT0 = fParam->GetMeanT0();
d0bcd1fb	284
bce12dc5	285
	286	for (Int_t in=0; in<12; in++)
	287	{
	288	timeCFD[in] = allData[in+1][0] ;
	289	timeCFD[in+12] = allData[in+56+1][0] ;
	290	timeLED[in] = allData[in+12+1][0] ;
	291	timeLED[in+12] = allData[in+68+1][0] ;
d0bcd1fb	292	AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
	293	in, timeCFD[in],timeCFD[in+12],timeLED[in],
	294	timeLED[in+12]));
bce12dc5	295	}
	296
	297	for (Int_t in=0; in<12; in++)
	298	{
d0bcd1fb	299	chargeQT0[in]=allData[2*in+25][0];
d0bcd1fb	300	chargeQT1[in]=allData[2*in+26][0];
bce12dc5	301	}
	302
	303	for (Int_t in=12; in<24; in++)
	304	{
d0bcd1fb	305	chargeQT0[in]=allData[2*in+57][0];
d0bcd1fb	306	chargeQT1[in]=allData[2*in+58][0];
bce12dc5	307	}
bce12dc5	308
d0bcd1fb	309	// } //cosmic with physics event
2e6a5ee0	310	for (Int_t in=0; in<24; in++)
c883fdf2	311	AliDebug(10, Form(" readed Raw %i %i %i %i %i",
c883fdf2	312	in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
776de217	313	onlineMean = allData[49][0];
2e6a5ee0	314
12e9daf9	315	Float_t time[24], adc[24];
2e6a5ee0	316	for (Int_t ipmt=0; ipmt<24; ipmt++) {
2e6a5ee0	317	if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
bce12dc5	318	//for simulated data
d0bcd1fb	319	//for physics data
	320	if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0)
	321	adc[ipmt] = chargeQT1[ipmt] - chargeQT0[ipmt];
	322	else
	323	adc[ipmt] = 0;
bce12dc5	324
539b9cb9	325
c555f418	326	time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
bce12dc5	327
d0bcd1fb	328	Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
539b9cb9	329	// time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt] ) ;
d0bcd1fb	330	AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
	331	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
	332	Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	333	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
	334	AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %i\n ",ampMip,qtMip, adc[ipmt]));
bce12dc5	335
d0bcd1fb	336	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	337	frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam
	338	frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
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	}
	362	if(besttimeA !=9999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
	363	if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
539b9cb9	364	AliDebug(5,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f #channel",
d0bcd1fb	365	pmtBestA,besttimeA, pmtBestC, besttimeC));
776de217	366	if(besttimeA <9999999 && besttimeC < 9999999 ){
9b83615d	367	timeDiff = ( besttimeA - besttimeC) *channelWidth;
d0bcd1fb	368	meanTime = Float_t((besttimeA + besttimeC)/2.);
776de217	369	timeclock = Float_t(meanTime - ref);
d0bcd1fb	370	vertex = meanVertex - c*(timeDiff)/2.; //+ (fdZonA - fdZonC)/2;
2e6a5ee0	371	}
776de217	372	} //if phys event
	373	frecpoints->SetT0clock(timeclock);
	374	frecpoints->SetVertex(vertex);
	375	frecpoints->SetMeanTime(Int_t(meanTime));
	376	frecpoints->SetOnlineMean(Int_t(onlineMean));
	377	AliDebug(5,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm meanVertex %f online mean %i ",timeDiff, meanTime,vertex,meanVertex, onlineMean));
	378
	379
bce12dc5	380	} // if (else )raw data
	381	recTree->Fill();
	382	if(frecpoints) delete frecpoints;
c41ceaac	383	}
bce12dc5	384
bce12dc5	385
c41ceaac	386	//____________________________________________________________
c41ceaac	387
d76c31f4	388	void AliT0Reconstructor::FillESD(TTree /digitsTree/, TTree clustersTree, AliESDEvent *pESD) const
dc7ca31d	389	{
	390
	391	/***************************************************
	392	Resonstruct digits to vertex position
	393	****************************************************/
	394
dc7ca31d	395	AliDebug(1,Form("Start FillESD T0"));
dc7ca31d	396
d76c31f4	397	TTree *treeR = clustersTree;
dc7ca31d	398
	399	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	400	if (!frecpoints) {
	401	AliError("Reconstruct Fill ESD >> no recpoints found");
	402	return;
	403	}
	404
	405	AliDebug(1,Form("Start FillESD T0"));
	406	TBranch *brRec = treeR->GetBranch("T0");
	407	if (brRec) {
	408	brRec->SetAddress(&frecpoints);
	409	}else{
f16935f7	410	AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d	411	return;
	412	}
	413
	414	brRec->GetEntry(0);
f16935f7	415	Float_t amp[24], time[24];
	416	Float_t zPosition = frecpoints -> GetVertex();
	417	Float_t timeStart = frecpoints -> GetMeanTime() ;
776de217	418	Float_t timeClock = frecpoints -> GetT0clock() ;
776de217	419	for ( Int_t i=0; i<24; i++) {
94c27e4f	420	time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
dc7ca31d	421	amp[i] = frecpoints -> GetAmp(i);
dc7ca31d	422	}
f16935f7	423	pESD->SetT0zVertex(zPosition); //vertex Z position
dc7ca31d	424	pESD->SetT0(timeStart); // interaction time
776de217	425	pESD->SetT0clock(timeClock); // interaction time with ref.point(spectrum)
dc7ca31d	426	pESD->SetT0time(time); // best TOF on each PMT
dc7ca31d	427	pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT
5325480c	428
f16935f7	429	AliDebug(1,Form(" Z position %f cm, T0 %f ps",zPosition , timeStart));
dc7ca31d	430
dc7ca31d	431	} // vertex in 3 sigma
	432
	433
	434
	435
	436
	437