[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;
    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 = ( besttimeC - besttimeA) *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 ){
12e9daf9	204	timeDiff = (besttimeC - besttimeA)*channelWidth;
d0bcd1fb	205	meanTime = Float_t((besttimeA + besttimeC)/2);// * channelWidth);
12e9daf9	206	// meanTime = (meanT0 - (besttimeA + besttimeC)/2) * channelWidth;
d0bcd1fb	207	vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
dc7ca31d	208	frecpoints->SetVertex(vertex);
dc7ca31d	209	frecpoints->SetMeanTime(Int_t(meanTime));
446d6ec4	210	//online mean
d0bcd1fb	211	frecpoints->SetOnlineMean(Int_t(onlineMean));
539b9cb9	212	AliDebug(10,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm online mean %i ps",timeDiff, meanTime,vertex, Int_t(onlineMean)));
dc7ca31d	213
dc7ca31d	214	}
b95e8d87	215
dc7ca31d	216	clustersTree->Fill();
bd375212	217
	218	delete timeCFD;
	219	delete timeLED;
	220	delete chargeQT0;
	221	delete chargeQT1;
dc7ca31d	222	}
	223
	224
c41ceaac	225	//_______________________________________________________________________
	226
	227	void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
	228	{
94c27e4f	229	// T0 raw ->
539b9cb9	230	//
	231	// reference amplitude and time ref. point from reco param
	232
	233	Int_t refAmp = GetRecoParam()->GetRefAmp();
	234	Int_t refPoint = GetRecoParam()->GetRefPoint();
	235
	236
e8ed1cd0	237	Int_t allData[110][5];
2e6a5ee0	238
e8ed1cd0	239	Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
2e6a5ee0	240
bce12dc5	241	for (Int_t i0=0; i0<105; i0++)
	242	{
	243	for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
	244	}
2e6a5ee0	245
bce12dc5	246	Float_t besttimeA=9999999;
	247	Float_t besttimeC=9999999;
	248	Int_t pmtBestA=99999;
	249	Int_t pmtBestC=99999;
	250	Float_t timeDiff=9999999, meanTime=0;
b95e8d87	251	Float_t meanVertex = fParam->GetMeanVertex();
29a60970	252
bce12dc5	253	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	254
	255	recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
2e6a5ee0	256
bce12dc5	257	AliDebug(10," before read data ");
	258	AliT0RawReader myrawreader(rawReader);
	259	if (!myrawreader.Next())
	260	AliDebug(1,Form(" no raw data found!!"));
	261	else
	262	{
	263	for (Int_t i=0; i<105; i++) {
	264	for (Int_t iHit=0; iHit<5; iHit++)
	265	{
	266	allData[i][iHit] = myrawreader.GetData(i,iHit);
	267	}
	268	}
	269
	270	Float_t channelWidth = fParam->GetChannelWidth() ;
	271
	272	// Int_t meanT0 = fParam->GetMeanT0();
d0bcd1fb	273
bce12dc5	274
	275	for (Int_t in=0; in<12; in++)
	276	{
	277	timeCFD[in] = allData[in+1][0] ;
	278	timeCFD[in+12] = allData[in+56+1][0] ;
	279	timeLED[in] = allData[in+12+1][0] ;
	280	timeLED[in+12] = allData[in+68+1][0] ;
d0bcd1fb	281	AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
	282	in, timeCFD[in],timeCFD[in+12],timeLED[in],
	283	timeLED[in+12]));
bce12dc5	284	}
	285
	286	for (Int_t in=0; in<12; in++)
	287	{
d0bcd1fb	288	chargeQT0[in]=allData[2*in+25][0];
d0bcd1fb	289	chargeQT1[in]=allData[2*in+26][0];
bce12dc5	290	}
	291
	292	for (Int_t in=12; in<24; in++)
	293	{
d0bcd1fb	294	chargeQT0[in]=allData[2*in+57][0];
d0bcd1fb	295	chargeQT1[in]=allData[2*in+58][0];
bce12dc5	296	}
bce12dc5	297
d0bcd1fb	298	// } //cosmic with physics event
2e6a5ee0	299	for (Int_t in=0; in<24; in++)
c883fdf2	300	AliDebug(10, Form(" readed Raw %i %i %i %i %i",
c883fdf2	301	in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
446d6ec4	302	Int_t onlineMean = allData[49][0];
2e6a5ee0	303
12e9daf9	304	Float_t time[24], adc[24];
2e6a5ee0	305	for (Int_t ipmt=0; ipmt<24; ipmt++) {
2e6a5ee0	306	if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
bce12dc5	307	//for simulated data
d0bcd1fb	308	//for physics data
	309	if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0)
	310	adc[ipmt] = chargeQT1[ipmt] - chargeQT0[ipmt];
	311	else
	312	adc[ipmt] = 0;
bce12dc5	313
539b9cb9	314
539b9cb9	315	time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
bce12dc5	316
d0bcd1fb	317	Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
539b9cb9	318	// time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt] ) ;
d0bcd1fb	319	AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
	320	ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
	321	Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	322	Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
	323	AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %i\n ",ampMip,qtMip, adc[ipmt]));
bce12dc5	324
d0bcd1fb	325	frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	326	frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam
	327	frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
bce12dc5	328
2e6a5ee0	329	}
	330	else {
	331	time[ipmt] = 0;
	332	adc[ipmt] = 0;
	333	}
	334	}
	335
	336	for (Int_t ipmt=0; ipmt<12; ipmt++){
	337	if(time[ipmt] > 1 ) {
	338	if(time[ipmt]<besttimeC){
	339	besttimeC=time[ipmt]; //timeC
	340	pmtBestC=ipmt;
	341	}
	342	}
	343	}
	344	for ( Int_t ipmt=12; ipmt<24; ipmt++){
	345	if(time[ipmt] > 1) {
	346	if(time[ipmt]<besttimeA) {
	347	besttimeA=time[ipmt]; //timeA
	348	pmtBestA=ipmt;}
	349	}
	350	}
	351	if(besttimeA !=9999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
	352	if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
539b9cb9	353	AliDebug(5,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f #channel",
d0bcd1fb	354	pmtBestA,besttimeA, pmtBestC, besttimeC));
d0bcd1fb	355	Float_t c = 0.0299792458; // cm/ps
2e6a5ee0	356	Float_t vertex = 99999;
2e6a5ee0	357	if(besttimeA <9999999 && besttimeC < 9999999 ){
12e9daf9	358	timeDiff = ( besttimeC - besttimeA) *channelWidth;
d0bcd1fb	359	meanTime = Float_t((besttimeA + besttimeC)/2.);
	360	onlineMean = onlineMean ;
	361	vertex = meanVertex - c*(timeDiff)/2.; //+ (fdZonA - fdZonC)/2;
2e6a5ee0	362	frecpoints->SetVertex(vertex);
2e6a5ee0	363	frecpoints->SetMeanTime(Int_t(meanTime));
1d7681da	364	frecpoints->SetOnlineMean(Int_t(onlineMean));
539b9cb9	365	AliDebug(5,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm meanVertex %f online mean %i ",timeDiff, meanTime,vertex,meanVertex, onlineMean));
2e6a5ee0	366
2e6a5ee0	367	}
bce12dc5	368	} // if (else )raw data
	369	recTree->Fill();
	370	if(frecpoints) delete frecpoints;
c41ceaac	371	}
bce12dc5	372
bce12dc5	373
c41ceaac	374	//____________________________________________________________
c41ceaac	375
d76c31f4	376	void AliT0Reconstructor::FillESD(TTree /digitsTree/, TTree clustersTree, AliESDEvent *pESD) const
dc7ca31d	377	{
	378
	379	/***************************************************
	380	Resonstruct digits to vertex position
	381	****************************************************/
	382
dc7ca31d	383	AliDebug(1,Form("Start FillESD T0"));
dc7ca31d	384
d76c31f4	385	TTree *treeR = clustersTree;
dc7ca31d	386
	387	AliT0RecPoint* frecpoints= new AliT0RecPoint ();
	388	if (!frecpoints) {
	389	AliError("Reconstruct Fill ESD >> no recpoints found");
	390	return;
	391	}
	392
	393	AliDebug(1,Form("Start FillESD T0"));
	394	TBranch *brRec = treeR->GetBranch("T0");
	395	if (brRec) {
	396	brRec->SetAddress(&frecpoints);
	397	}else{
f16935f7	398	AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d	399	return;
	400	}
	401
	402	brRec->GetEntry(0);
f16935f7	403	Float_t amp[24], time[24];
	404	Float_t zPosition = frecpoints -> GetVertex();
	405	Float_t timeStart = frecpoints -> GetMeanTime() ;
	406	for ( Int_t i=0; i<24; i++) {
94c27e4f	407	time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
dc7ca31d	408	amp[i] = frecpoints -> GetAmp(i);
dc7ca31d	409	}
f16935f7	410	pESD->SetT0zVertex(zPosition); //vertex Z position
dc7ca31d	411	pESD->SetT0(timeStart); // interaction time
	412	pESD->SetT0time(time); // best TOF on each PMT
	413	pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT
5325480c	414
f16935f7	415	AliDebug(1,Form(" Z position %f cm, T0 %f ps",zPosition , timeStart));
dc7ca31d	416
dc7ca31d	417	} // vertex in 3 sigma
	418
	419
	420
	421
	422
	423