[u/mrichter/AliRoot.git] / T0 / AliT0Reconstructor.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */
/*********************************************************************
 *  T0 reconstruction and filling ESD
 *  - reconstruct mean time (interation time) 
 *  - vertex position
 *  -  multiplicity
 ********************************************************************/

#include <AliESDEvent.h>
#include "AliLog.h"
#include "AliT0RecPoint.h"
#include "AliRawReader.h"
#include "AliT0RawReader.h"
#include "AliT0digit.h"
#include "AliT0Reconstructor.h"
#include "AliT0Parameters.h"
#include "AliT0Calibrator.h"

#include <TArrayI.h>
#include <TGraph.h>
#include <TMath.h>
#include <Riostream.h>

ClassImp(AliT0Reconstructor)

  AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
					     fdZonA(0),
					     fdZonC(0),
					     fZposition(0),
					     fParam(NULL),
					     fAmpLEDrec(),
					     fQTC(0),
					     fAmpLED(0),
					     fCalib()
{
  //constructor

  fParam = AliT0Parameters::Instance();
  fParam->Init();
 
  for (Int_t i=0; i<24; i++){
        TGraph* gr = fParam ->GetAmpLEDRec(i);
	if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ; 
	  TGraph* gr1 = fParam ->GetAmpLED(i);
	  if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ; 
	  TGraph* gr2 = fParam ->GetQTC(i);
	  if (gr2) fQTC.AddAtAndExpand(gr2,i) ; 	
  }

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


  fCalib = new AliT0Calibrator();

}

//_____________________________________________________________________________
void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
  
{
  // T0 digits reconstruction
  Int_t refAmp = GetRecoParam()->GetRefAmp();
  
  TArrayI * timeCFD = new TArrayI(24); 
  TArrayI * timeLED = new TArrayI(24); 
  TArrayI * chargeQT0 = new TArrayI(24); 
  TArrayI * chargeQT1 = new TArrayI(24); 

 
  Float_t channelWidth = fParam->GetChannelWidth() ;  
  Float_t meanVertex = fParam->GetMeanVertex();
  Float_t c = 0.0299792; // cm/ps
  Double32_t vertex = 9999999;
  Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;

  
  AliDebug(1,Form("Start DIGITS reconstruction "));
  

  TBranch *brDigits=digitsTree->GetBranch("T0");
  AliT0digit *fDigits = new AliT0digit() ;
  if (brDigits) {
    brDigits->SetAddress(&fDigits);
  }else{
    AliError(Form("EXEC Branch T0 digits not found"));
     return;
  }
  
  digitsTree->GetEvent(0);
  digitsTree->GetEntry(0);
  brDigits->GetEntry(0);
  fDigits->GetTimeCFD(*timeCFD);
  fDigits->GetTimeLED(*timeLED);
  fDigits->GetQT0(*chargeQT0);
  fDigits->GetQT1(*chargeQT1);
  Int_t onlineMean =  fDigits->MeanTime();

  Bool_t tr[5];
  for (Int_t i=0; i<5; i++) tr[i]=false; 
  
  Double32_t besttimeA=999999;
  Double32_t besttimeC=999999;
  Int_t pmtBestA=99999;
  Int_t pmtBestC=99999;
  
  AliT0RecPoint* frecpoints= new AliT0RecPoint ();
  clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
  
  Float_t time[24], adc[24];
  for (Int_t ipmt=0; ipmt<24; ipmt++) {
    if(timeCFD->At(ipmt)>0 ){
     if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)  
	adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
      else
	adc[ipmt] = 0;
      
     // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt],  timeCFD->At(ipmt)) ;
	     
      Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
      time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl),  timeCFD->At(ipmt) ) ;
      AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
		       ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));

      Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
      Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
      AliDebug(10,Form("  Amlitude in MIPS LED %f ,  QTC %f in channels %i\n ",ampMip,qtMip, adc[ipmt]));
      
      frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
      frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam 
      frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
      
    }
    else {
      time[ipmt] = 0;
      adc[ipmt] = 0;
    }
  }
  
  for (Int_t ipmt=0; ipmt<12; ipmt++){
    if(time[ipmt] > 1 ) {
      if(time[ipmt]<besttimeC){
	besttimeC=time[ipmt]; //timeC
	pmtBestC=ipmt;
      }
    }
  }
  for ( Int_t ipmt=12; ipmt<24; ipmt++){
    if(time[ipmt] > 1) {
      if(time[ipmt]<besttimeA) {
	besttimeA=time[ipmt]; //timeA
        pmtBestA=ipmt;}
    }
  }
  if(besttimeA < 999999) {
    frecpoints->SetTimeBestA(Int_t(besttimeA *channelWidth));
    tr[1]=true;
  }
  if( besttimeC < 999999 ) {
    frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth));
    tr[2]=true;
  }
  AliDebug(10,Form(" besttimeA %f ch,  besttimeC %f ch",besttimeA, besttimeC));
  if(besttimeA <999999 && besttimeC < 999999 ){
    //    timeDiff = (besttimeC - besttimeA)*channelWidth;
    timeDiff = (besttimeA - besttimeC)*channelWidth;
    meanTime = (besttimeA + besttimeC)/2;// * channelWidth); 
    timeclock = meanTime *channelWidth ;
    vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
    tr[0]=true; 
  }
  frecpoints->SetVertex(vertex);
  frecpoints->SetMeanTime(meanTime);
  frecpoints->SetT0clock(timeclock);
  frecpoints->SetT0Trig(tr);

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

  //online mean
  frecpoints->SetOnlineMean(Int_t(onlineMean));
  AliDebug(10,Form("  timeDiff %i #channel,  meanTime %i #channel, vertex %f cm online mean %i timeclock %i ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
  
  
  clustersTree->Fill();

  delete timeCFD;
  delete timeLED;
  delete chargeQT0; 
  delete chargeQT1; 
}


//_______________________________________________________________________

void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
{
  // T0 raw ->
  //
  // reference amplitude and time ref. point from reco param

  Int_t refAmp = GetRecoParam()->GetRefAmp();
  Int_t refPoint = GetRecoParam()->GetRefPoint();
  Float_t latencyL1 = GetRecoParam()->GetLatencyL1();
  Float_t latencyHPTDC = GetRecoParam()->GetLatencyHPTDC();
  Float_t vertexshift = GetRecoParam()->GetVertexShift();

  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;

  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;
  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=0;
      if (refPoint>0) 
      ref = allData[refPoint][0]-5000;

      Float_t channelWidth = fParam->GetChannelWidth() ;  
      
      //       Int_t meanT0 = fParam->GetMeanT0();
       
	  
	  for (Int_t in=0; in<12; in++)  
	    {
	      timeCFD[in] = allData[in+1][0] ;
	      timeCFD[in+12] = allData[in+56+1][0] ;
	      timeLED[in] = allData[in+12+1][0] ;
	      timeLED[in+12] = allData[in+68+1][0] ;
	      AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
				in, timeCFD[in],timeCFD[in+12],timeLED[in], 
				timeLED[in+12]));   
	    }
	  
	  for (Int_t in=0; in<12;  in++)
	    {
	      chargeQT0[in]=allData[2*in+25][0];
	      chargeQT1[in]=allData[2*in+26][0];
	    }
	  
	   for (Int_t in=12; in<24;  in++)
	     {
	       chargeQT0[in]=allData[2*in+57][0];
	       chargeQT1[in]=allData[2*in+58][0];
	     }
	   
	   //	 } //cosmic with physics event
       for (Int_t in=0; in<24; in++)  
	 AliDebug(10, Form(" readed Raw %i %i %i %i %i",
			   in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
        onlineMean = allData[49][0];       
       
       Double32_t time[24], adc[24];
       for (Int_t ipmt=0; ipmt<24; ipmt++) {
	 if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
	   //for simulated data
	     //for physics  data
	   if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0)  
	     adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
	   else
	     adc[ipmt] = 0;
	   

	   //	   time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
	   
      	   Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
	   time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
	   AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
			    ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
	   Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
	   Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
	   AliDebug(10,Form("  Amlitude in MIPS LED %f ; QTC %f;  in channels %i\n ",ampMip,qtMip, adc[ipmt]));
	     
	   //	   frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
	   frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
	   frecpoints->SetAmpLED(ipmt, Double32_t( qtMip)); //for cosmic &pp beam 
	   frecpoints->SetAmp(ipmt, Double32_t(ampMip));
	     
	 }
	 else {
	   time[ipmt] = 0;
	   adc[ipmt] = 0;
	 }
       }
       
       for (Int_t ipmt=0; ipmt<12; ipmt++){
	 if(time[ipmt] > 1 ) {
	   if(time[ipmt]<besttimeC){
	     besttimeC=time[ipmt]; //timeC
	     pmtBestC=ipmt;
	   }
	 }
       }
       for ( Int_t ipmt=12; ipmt<24; ipmt++){
	 if(time[ipmt] > 1) {
	   if(time[ipmt]<besttimeA) {
	     besttimeA=time[ipmt]; //timeA
	     pmtBestA=ipmt;}
	 }
       }
       if(besttimeA < 999999) 
	 frecpoints->SetTimeBestA( 0.001 * besttimeA * channelWidth - latencyHPTDC + latencyL1);
       if( besttimeC < 999999 ) 
	 frecpoints->SetTimeBestC(0.001 * besttimeC * channelWidth - latencyHPTDC + latencyL1);
       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;
	 timeclock =  0.001 * channelWidth * Float_t(besttimeA+besttimeC)/2. - latencyHPTDC + latencyL1;  
	 meanTime = (besttimeA+besttimeC-2.*Float_t(ref))/2.;
	 vertex =  meanVertex -  c*(timeDiff)/2. + vertexshift; //+ (fdZonA - fdZonC)/2; 
	}
      }  //if phys event       
      AliDebug(5,Form("  timeDiff %f #channel,  meanTime %f #channel, TOFmean%f  vertex %f cm meanVertex %f online mean %i \n",timeDiff, meanTime,timeclock, vertex,meanVertex, onlineMean));
      frecpoints->SetT0clock(timeclock);
      frecpoints->SetVertex(vertex);
      frecpoints->SetMeanTime(meanTime);
      frecpoints->SetOnlineMean(Int_t(onlineMean));
	// Set triggers
      
      Bool_t tr[5];
      Int_t trchan[5]= {50,51,52,55,56};
      for (Int_t i=0; i<5; i++) tr[i]=false; 
      for (Int_t itr=0; itr<5; itr++) {
	if(allData[trchan[itr]][0]>0) tr[itr]=true;
	frecpoints->SetT0Trig(tr);
      }
      for (Int_t i=0; i<5; i++) 
	printf(" T0 trigers %i ",tr[i]);
      printf(" \n ");
      
      
      } // if (else )raw data
      recTree->Fill();
      if(frecpoints) delete frecpoints;
    }
  
  
  //____________________________________________________________
  
  void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
  {

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