[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();
  Int_t refPoint = GetRecoParam()->GetRefPoint();
 
  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();
  
  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();

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