[u/mrichter/AliRoot.git] / PHOS / AliPHOSDApi0mip.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$ */

// --
// --
// Implementation for TTree output in PHOS DA
// for calibrating energy by pi0 and MIP.
// --
// -- Author: Hisayuki Torii (Hiroshima Univ.)
// --

#include <stdio.h>
#include <iostream>
#include <math.h>
#include "AliLog.h"
#include "AliPHOSDApi0mip.h"
ClassImp(AliPHOSDApi0mip)
//----------------------------------------------------------------
AliPHOSDApi0mip::AliPHOSDApi0mip(int module,int iterid,const char* fopt) :
  TNamed(), fCreateTree(false), fCreateHist(false), fMod(0), fIterId(0),
  fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
  fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
  fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
  // Constructor

  char hname[1024], htitle[1024];

  TString shname="AliPHOSDApi0mip_mod%d_iter%d";
  snprintf(hname,shname.Length(),shname.Data(),module,iterid);
  SetName(hname);

  TString shtitle="PHOS MIP/pi0 Calibration DA for Module:%d Iteration:%d";
  snprintf(htitle,shtitle.Length(),shtitle.Data(),module,iterid);
  SetTitle(htitle);

  fMod = module;
  fIterId = iterid;
  fEvent = 0;
  fEventClustered = false;
  fTime = 0;
  fCreateTree = false;
  fCreateHist = false;

  char fname[1024];

  TString sfname="AliPHOSDApi0mip_mod%d.root";
  snprintf(fname,sfname.Length(),sfname.Data(),module);
  fTFile = TFile::Open(fname,fopt);
}
//----------------------------------------------------------------
AliPHOSDApi0mip::AliPHOSDApi0mip(const AliPHOSDApi0mip& da):
  TNamed(da), fCreateTree(false), fCreateHist(false), fMod(da.fMod), fIterId(da.fIterId),
  fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
  fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
  fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
  // Copy Constructor

  char fname[1024], hname[1024], htitle[1024];;

  TString sfname="%s.root";
  snprintf(fname,sfname.Length(),sfname.Data(),GetName());
  fTFile = TFile::Open(fname,"RECREATE");
  fEvent = 0;
  fEventClustered = false;
  fTime = 0;

  if( da.fCreateHist ){
    fCreateHist = true;
    fH1Time = (TH1I*) da.fH1Time->Clone();
    fH1DigitNum = (TH1F*) da.fH1DigitNum->Clone();
    fH1ClusterNum = (TH1F*) da.fH1ClusterNum->Clone();
    fH2EneDigitId = (TH2F*) da.fH2EneDigitId->Clone();
    fH2MipDigitId = (TH2F*) da.fH2MipDigitId->Clone();
    fH2Pi0DigitId = (TH2F*) da.fH2Pi0DigitId->Clone();
    fH3Pi0AsymPt = (TH3F*) da.fH3Pi0AsymPt->Clone();
  } else {
    fCreateHist = false;
    fH1Time = 0;
    fH1DigitNum = 0;
    fH1ClusterNum = 0;
    fH2EneDigitId = 0;
    fH2MipDigitId = 0;
    fH2Pi0DigitId = 0;
    fH3Pi0AsymPt = 0;
  }
  if( da.fCreateTree ){
    // Create new ttree instead of copy.
    fCreateTree = true;

    TString shname="tevt_mod%d_iter%d";
    snprintf(hname,shname.Length(),shname.Data(),fMod,fIterId);

    TString shtitle="Calibration for Module:%d Iteration:%d";
    snprintf(htitle,shtitle.Length(),shtitle.Data(),fMod,fIterId);

    fTTree = new TTree(hname,htitle);
    fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
  } else {
    fCreateTree = false;
    fTTree = 0;
  }
}
//-------------------------------------------------------------------
AliPHOSDApi0mip& AliPHOSDApi0mip::operator= (const AliPHOSDApi0mip&)
{
  // Operator= is not implemented yet
  AliFatal("Operator= is not implemented");  
  return *this;
}
//-------------------------------------------------------------------
AliPHOSDApi0mip::~AliPHOSDApi0mip(){
  // Destructor

  //std::cout<<" DEBUGDEBUG :: AliPHOSDApi0mip:: deconstructor "<<std::endl;
  if( fCreateHist ){
    fH1Time->Write(0,TObject::kOverwrite);
    fH1DigitNum->Write(0,TObject::kOverwrite);
    fH1ClusterNum->Write(0,TObject::kOverwrite);
    fH2EneDigitId->Write(0,TObject::kOverwrite);
    fH2MipDigitId->Write(0,TObject::kOverwrite);
    fH2Pi0DigitId->Write(0,TObject::kOverwrite);
    fH3Pi0AsymPt->Write(0,TObject::kOverwrite);
    delete fH1Time;
    delete fH1DigitNum;
    delete fH1ClusterNum;
    delete fH2EneDigitId;
    delete fH2MipDigitId;
    delete fH2Pi0DigitId;
    delete fH3Pi0AsymPt;
  }
  if( fCreateTree ){
    fTTree->Write();
  }
  if( fEvent ) delete fEvent;
  fTFile->Close();
  delete fTFile;
}
//-------------------------------------------------------------------
void AliPHOSDApi0mip::FillDigit(float adc,int row,int col){
  // Put new digit information in event
  //

  if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
  float energy = adc * 0.005; // 5MeV/ch is fixed.
  if( energy > 0.020 ){  // Threshold at 20MeV.
    //std::cout<<" AliPHOSDApi0mip::AppendDigit() DEBUG:: (row,col)=("<<row<<","<<col<<") adc="<<adc<<" energy="<<energy<<std::endl;
    fEvent->Fill(energy,row,col);
  }
}
//-------------------------------------------------------------------
void AliPHOSDApi0mip::NewEvent(){
  // Initiator must be called in the beginning of event
  //

  if(fEvent) fEvent->Reset();
  fEventClustered = false;
}
//----------------------------------------------------------------
bool AliPHOSDApi0mip::CreateTree(){
  if(! fTFile ) {
    return false;
  }
  fTFile->cd();
  char hname[100], htitle[100];
  snprintf(hname,100,"trevt_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"Calibration for Module:%d Iteration:%d",fMod,fIterId);
  fTTree = new TTree(hname,htitle);
  fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
  fCreateTree = true;
  return true;
}
//----------------------------------------------------------------
bool AliPHOSDApi0mip::CreateHist(){
  // Create histogram routine called by the constructor only.

  char hname[100], htitle[100];

  snprintf(hname,100,"h1_time_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"Time : Mod:%d Iter:%d",fMod,fIterId);
  fH1Time = (TH1I*) gDirectory->Get(hname);
  if( fH1Time>0 ){
    std::cout<<" AliPHOSDApi0mip:Warning!! Output object already exist : "<<fH1Time->GetName()<<std::endl;
    return false;
  }
  fH1Time = new TH1I(hname,htitle,2,0,2);

  snprintf(hname,100,"h1_digitnum_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"Number of Digits : Mod:%d Iter:%d",fMod,fIterId);
  fH1DigitNum = new TH1F(hname,htitle,100,0,100);

  snprintf(hname,100,"h1_clusternum_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"Number of Clusters : Mod:%d Iter:%d",fMod,fIterId);
  fH1ClusterNum = new TH1F(hname,htitle,100,0,100);

  snprintf(hname,100,"h2_pi0digitid_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"PHOS pi0 mass vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
  fH2Pi0DigitId = new TH2F(hname,htitle,17920,0,17920,120,0,0.3);
      
  snprintf(hname,100,"h2_mipdigitid_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
  fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
  //fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,50,0.5,1.5);

  snprintf(hname,100,"h2_enedigitid_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
  //fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
  fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,50,0,5.0);

  snprintf(hname,100,"h3_pi0asympt_mod%d_iter%d",fMod,fIterId);
  snprintf(htitle,100,"PHOS pi0 mass vs pT vs Asym : Mod:%d Iter:%d",fMod,fIterId);
  fH3Pi0AsymPt = new TH3F(hname,htitle,20,0,1,20,0,10,200,0,1);

  fCreateHist = true;
  return true;
}
//-------------------------------------------------------------------
void AliPHOSDApi0mip::FillTree(AliPHOSDATreeEvent* event){
  // Fill all information (AliPHOSDATreeEvent) into TTree
  // if event==NULL, the obtained information (=fEvent) is used.

  if(! fCreateTree ) CreateTree();
  if( event == 0 ){
    if( fEvent == 0 ){
      fEvent = new AliPHOSDATreeEvent();
    }
    // No need to make clustering for tree saving.
    //if(! fEventClustered ){
    //fEvent->ExecuteClustering();
    //fEventClustered = true;
    //}
    fEvent->SetTime(fTime);
    fTTree->Fill();
  } else {
    if( fEvent ) delete fEvent;
    fEvent = event;
    fTTree->Fill();
    fEvent = 0;
  }
}
//-------------------------------------------------------------------
void AliPHOSDApi0mip::FillHist(AliPHOSDATreeEvent* event){
  // Fill all information (AliPHOSDATreeEvent) into histograms
  // if event==NULL, the obtained information (=fEvent) is used.

  if(! fCreateHist ) CreateHist();
  if( event == 0 ){
    if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
    if(! fEventClustered ){
      fEvent->ExecuteClustering();
      fEventClustered = true;
    }
    fEvent->SetTime(fTime);
    event = fEvent;
  }

  // -- Constant
  int nclt1, nclt2;
  int ndigits;
  int ndigitsall = 0;
  AliPHOSDATreeDigit digit1, digit2;
  const float kDistanceFromVertex = 460.;
  float dist1, dist2, px, py;
  float mass, cosine, pt, asym;

  // -- Filling into Time Histogram
  double start = fH1Time->GetBinContent(1);
  double stop = fH1Time->GetBinContent(2);
  double currenttime = (double) (event->GetTime());
  if( start == 0 && stop == 0 ) start = stop = currenttime;
  if( currenttime < start ) start = currenttime;
  if( currenttime > stop ) stop = currenttime;
  fH1Time->SetBinContent(1,start);
  fH1Time->SetBinContent(2,stop);

  // -- Start looping clusters
  nclt1 = event->GetNClusters();
  fH1ClusterNum->Fill(nclt1);
  //std::cout<<" AliPHOSDApi0mip::FillHisto() DEBUG:: Number of clusters: "<<nclt1<<std::endl;
  while( nclt1-- ){
    AliPHOSDATreeCluster& clt1 = event->GetCluster(nclt1);
    ndigits = clt1.GetNDigits();
    while( ndigits-- ){
      AliPHOSDATreeDigit& digit = clt1.GetDigit(ndigits);
      fH2EneDigitId->Fill(digit.GetDigitId(),digit.GetEnergy());
      ndigitsall++;
    }
    digit1 = clt1.GetMaxDigit();
    if( clt1.GetEnergy() > 0 && clt1.GetNDigits() >= 2 ){
      fH2MipDigitId->Fill(digit1.GetDigitId(),clt1.GetEnergy());
    }
    if( clt1.GetEnergy() > 0.2 ){
      nclt2 = nclt1;
      while( nclt2-- ){
	AliPHOSDATreeCluster& clt2 = event->GetCluster(nclt2);
	digit2 = clt2.GetMaxDigit();
	if( clt2.GetEnergy() > 0.2 ){
	  // 
	  // Following approximation is applied in order to avoid many calling "sqrt" function.
	  // sqrt( 1 + x^2 ) =~ 1 + 0.5x^2 (for x<<1)
	  //
	  //cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex)
	  //	    / ( kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*2.42/kDistanceFromVertex )
	  //	    / ( kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*2.42/kDistanceFromVertex );
	  //
	  

	  dist1 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*4.84 );
	  dist2 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 );
	  cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex)
	    / dist1 / dist2;
	  mass = sqrt( 2.0 * clt1.GetEnergy() * clt2.GetEnergy() * (1 - cosine) );
	  //std::cout<<" DEBUG: dist1:"<<dist1<<" dist2:"<<dist2<<" cosine:"<<cosine<<" c1:"<<clt1.GetEnergy()<<" c2:"<<clt2.GetEnergy()<<" mass:"<<mass<<std::endl;
	  //std::cout<<"      : clt1(row,col)=("<<clt1.GetRow()<<","<<clt1.GetCol()<<") clt2(row,col)=("<<clt2.GetRow()<<","<<clt2.GetCol()<<") "<<std::endl;
	  asym = fabs(clt1.GetEnergy() - clt2.GetEnergy())/(clt1.GetEnergy() + clt2.GetEnergy());
	  py = clt1.GetEnergy() * kDistanceFromVertex / dist1 + clt2.GetEnergy() * kDistanceFromVertex / dist2;
	  px = clt1.GetEnergy() * (clt1.GetCol()-27.5)*2.42/dist1 + clt2.GetEnergy() * (clt2.GetCol()-27.5)*2.42/ dist2;
	  pt = sqrt(px*px+py*py);
	  fH3Pi0AsymPt->Fill(asym,pt,mass);
	  if( digit1.IsValid() && digit1.IsValid() &&
	      digit1.GetDigitId() != digit2.GetDigitId() &&
	      //pt > 1.5 && clt1.GetEnergy() > 0.3 && clt2.GetEnergy() > 0.3 ){
	      //pt > 1.3 && clt1.GetEnergy() > 0.25 && clt2.GetEnergy() > 0.25 ){
	      pt > 1.0 && clt1.GetEnergy() > 0.2 && clt2.GetEnergy() > 0.2 ){
	    fH2Pi0DigitId->Fill(digit1.GetDigitId(),mass);
	    fH2Pi0DigitId->Fill(digit2.GetDigitId(),mass);
	  }
	} // End of energy cut
      } // End of second cluster loop (clt2)
    } // End of energy cut
  } // End of first cluster loop (clt1)
  fH1DigitNum->Fill(ndigitsall);
  //
}
//-------------------------------------------------------------------
void AliPHOSDApi0mip::Print(Option_t *option) const
{
  // Print Out

  //fTFile->ls();
  //fTTree->Print();
  if( fEvent ){
    //fEvent->ExecuteClustering();
    fEvent->Print(option);
  }
}
//-------------------------------------------------------------------
Commit	Line	Data
c335b3c5	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	**************************************************************************/
e7386040	15	/* $Id$ */
c335b3c5	16
12f6dd57	17	// --
	18	// --
	19	// Implementation for TTree output in PHOS DA
	20	// for calibrating energy by pi0 and MIP.
	21	// --
	22	// -- Author: Hisayuki Torii (Hiroshima Univ.)
	23	// --
	24
	25	#include <stdio.h>
	26	#include <iostream>
	27	#include <math.h>
c335b3c5	28	#include "AliLog.h"
12f6dd57	29	#include "AliPHOSDApi0mip.h"
	30	ClassImp(AliPHOSDApi0mip)
	31	//----------------------------------------------------------------
af695ce2	32	AliPHOSDApi0mip::AliPHOSDApi0mip(int module,int iterid,const char* fopt) :
12f6dd57	33	TNamed(), fCreateTree(false), fCreateHist(false), fMod(0), fIterId(0),
	34	fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
	35	fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
	36	fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
	37	// Constructor
	38
	39	char hname[1024], htitle[1024];
6ba1dd81	40
	41	TString shname="AliPHOSDApi0mip_mod%d_iter%d";
	42	snprintf(hname,shname.Length(),shname.Data(),module,iterid);
12f6dd57	43	SetName(hname);
6ba1dd81	44
	45	TString shtitle="PHOS MIP/pi0 Calibration DA for Module:%d Iteration:%d";
	46	snprintf(htitle,shtitle.Length(),shtitle.Data(),module,iterid);
12f6dd57	47	SetTitle(htitle);
	48
	49	fMod = module;
	50	fIterId = iterid;
	51	fEvent = 0;
	52	fEventClustered = false;
	53	fTime = 0;
	54	fCreateTree = false;
	55	fCreateHist = false;
	56
	57	char fname[1024];
6ba1dd81	58
	59	TString sfname="AliPHOSDApi0mip_mod%d.root";
	60	snprintf(fname,sfname.Length(),sfname.Data(),module);
12f6dd57	61	fTFile = TFile::Open(fname,fopt);
	62	}
	63	//----------------------------------------------------------------
	64	AliPHOSDApi0mip::AliPHOSDApi0mip(const AliPHOSDApi0mip& da):
	65	TNamed(da), fCreateTree(false), fCreateHist(false), fMod(da.fMod), fIterId(da.fIterId),
	66	fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
	67	fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
	68	fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
	69	// Copy Constructor
	70
	71	char fname[1024], hname[1024], htitle[1024];;
6ba1dd81	72
	73	TString sfname="%s.root";
	74	snprintf(fname,sfname.Length(),sfname.Data(),GetName());
12f6dd57	75	fTFile = TFile::Open(fname,"RECREATE");
	76	fEvent = 0;
	77	fEventClustered = false;
	78	fTime = 0;
	79
	80	if( da.fCreateHist ){
	81	fCreateHist = true;
	82	fH1Time = (TH1I*) da.fH1Time->Clone();
	83	fH1DigitNum = (TH1F*) da.fH1DigitNum->Clone();
	84	fH1ClusterNum = (TH1F*) da.fH1ClusterNum->Clone();
	85	fH2EneDigitId = (TH2F*) da.fH2EneDigitId->Clone();
	86	fH2MipDigitId = (TH2F*) da.fH2MipDigitId->Clone();
	87	fH2Pi0DigitId = (TH2F*) da.fH2Pi0DigitId->Clone();
	88	fH3Pi0AsymPt = (TH3F*) da.fH3Pi0AsymPt->Clone();
	89	} else {
	90	fCreateHist = false;
	91	fH1Time = 0;
	92	fH1DigitNum = 0;
	93	fH1ClusterNum = 0;
	94	fH2EneDigitId = 0;
	95	fH2MipDigitId = 0;
	96	fH2Pi0DigitId = 0;
	97	fH3Pi0AsymPt = 0;
	98	}
	99	if( da.fCreateTree ){
	100	// Create new ttree instead of copy.
	101	fCreateTree = true;
6ba1dd81	102
	103	TString shname="tevt_mod%d_iter%d";
	104	snprintf(hname,shname.Length(),shname.Data(),fMod,fIterId);
	105
	106	TString shtitle="Calibration for Module:%d Iteration:%d";
	107	snprintf(htitle,shtitle.Length(),shtitle.Data(),fMod,fIterId);
	108
12f6dd57	109	fTTree = new TTree(hname,htitle);
	110	fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
	111	} else {
	112	fCreateTree = false;
	113	fTTree = 0;
	114	}
	115	}
	116	//-------------------------------------------------------------------
c335b3c5	117	AliPHOSDApi0mip& AliPHOSDApi0mip::operator= (const AliPHOSDApi0mip&)
	118	{
	119	// Operator= is not implemented yet
	120	AliFatal("Operator= is not implemented");
	121	return *this;
	122	}
	123	//-------------------------------------------------------------------
12f6dd57	124	AliPHOSDApi0mip::~AliPHOSDApi0mip(){
	125	// Destructor
	126
	127	//std::cout<<" DEBUGDEBUG :: AliPHOSDApi0mip:: deconstructor "<<std::endl;
	128	if( fCreateHist ){
	129	fH1Time->Write(0,TObject::kOverwrite);
	130	fH1DigitNum->Write(0,TObject::kOverwrite);
	131	fH1ClusterNum->Write(0,TObject::kOverwrite);
	132	fH2EneDigitId->Write(0,TObject::kOverwrite);
	133	fH2MipDigitId->Write(0,TObject::kOverwrite);
	134	fH2Pi0DigitId->Write(0,TObject::kOverwrite);
	135	fH3Pi0AsymPt->Write(0,TObject::kOverwrite);
	136	delete fH1Time;
	137	delete fH1DigitNum;
	138	delete fH1ClusterNum;
	139	delete fH2EneDigitId;
	140	delete fH2MipDigitId;
	141	delete fH2Pi0DigitId;
	142	delete fH3Pi0AsymPt;
	143	}
	144	if( fCreateTree ){
	145	fTTree->Write();
	146	}
	147	if( fEvent ) delete fEvent;
	148	fTFile->Close();
	149	delete fTFile;
	150	}
	151	//-------------------------------------------------------------------
	152	void AliPHOSDApi0mip::FillDigit(float adc,int row,int col){
	153	// Put new digit information in event
	154	//
	155
	156	if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
	157	float energy = adc * 0.005; // 5MeV/ch is fixed.
	158	if( energy > 0.020 ){ // Threshold at 20MeV.
	159	//std::cout<<" AliPHOSDApi0mip::AppendDigit() DEBUG:: (row,col)=("<<row<<","<<col<<") adc="<<adc<<" energy="<<energy<<std::endl;
	160	fEvent->Fill(energy,row,col);
	161	}
	162	}
	163	//-------------------------------------------------------------------
	164	void AliPHOSDApi0mip::NewEvent(){
	165	// Initiator must be called in the beginning of event
	166	//
	167
	168	if(fEvent) fEvent->Reset();
	169	fEventClustered = false;
	170	}
	171	//----------------------------------------------------------------
	172	bool AliPHOSDApi0mip::CreateTree(){
	173	if(! fTFile ) {
	174	return false;
	175	}
	176	fTFile->cd();
14895588	177	char hname[100], htitle[100];
	178	snprintf(hname,100,"trevt_mod%d_iter%d",fMod,fIterId);
	179	snprintf(htitle,100,"Calibration for Module:%d Iteration:%d",fMod,fIterId);
12f6dd57	180	fTTree = new TTree(hname,htitle);
	181	fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
	182	fCreateTree = true;
	183	return true;
	184	}
	185	//----------------------------------------------------------------
	186	bool AliPHOSDApi0mip::CreateHist(){
	187	// Create histogram routine called by the constructor only.
	188
14895588	189	char hname[100], htitle[100];
12f6dd57	190
14895588	191	snprintf(hname,100,"h1_time_mod%d_iter%d",fMod,fIterId);
14895588	192	snprintf(htitle,100,"Time : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	193	fH1Time = (TH1I*) gDirectory->Get(hname);
	194	if( fH1Time>0 ){
	195	std::cout<<" AliPHOSDApi0mip:Warning!! Output object already exist : "<<fH1Time->GetName()<<std::endl;
	196	return false;
	197	}
	198	fH1Time = new TH1I(hname,htitle,2,0,2);
	199
14895588	200	snprintf(hname,100,"h1_digitnum_mod%d_iter%d",fMod,fIterId);
14895588	201	snprintf(htitle,100,"Number of Digits : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	202	fH1DigitNum = new TH1F(hname,htitle,100,0,100);
12f6dd57	203
14895588	204	snprintf(hname,100,"h1_clusternum_mod%d_iter%d",fMod,fIterId);
14895588	205	snprintf(htitle,100,"Number of Clusters : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	206	fH1ClusterNum = new TH1F(hname,htitle,100,0,100);
12f6dd57	207
14895588	208	snprintf(hname,100,"h2_pi0digitid_mod%d_iter%d",fMod,fIterId);
14895588	209	snprintf(htitle,100,"PHOS pi0 mass vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	210	fH2Pi0DigitId = new TH2F(hname,htitle,17920,0,17920,120,0,0.3);
12f6dd57	211
14895588	212	snprintf(hname,100,"h2_mipdigitid_mod%d_iter%d",fMod,fIterId);
14895588	213	snprintf(htitle,100,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	214	fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
	215	//fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,50,0.5,1.5);
	216
14895588	217	snprintf(hname,100,"h2_enedigitid_mod%d_iter%d",fMod,fIterId);
14895588	218	snprintf(htitle,100,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	219	//fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
	220	fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,50,0,5.0);
	221
14895588	222	snprintf(hname,100,"h3_pi0asympt_mod%d_iter%d",fMod,fIterId);
14895588	223	snprintf(htitle,100,"PHOS pi0 mass vs pT vs Asym : Mod:%d Iter:%d",fMod,fIterId);
12f6dd57	224	fH3Pi0AsymPt = new TH3F(hname,htitle,20,0,1,20,0,10,200,0,1);
	225
	226	fCreateHist = true;
	227	return true;
	228	}
	229	//-------------------------------------------------------------------
	230	void AliPHOSDApi0mip::FillTree(AliPHOSDATreeEvent* event){
	231	// Fill all information (AliPHOSDATreeEvent) into TTree
	232	// if event==NULL, the obtained information (=fEvent) is used.
	233
	234	if(! fCreateTree ) CreateTree();
	235	if( event == 0 ){
	236	if( fEvent == 0 ){
	237	fEvent = new AliPHOSDATreeEvent();
	238	}
	239	// No need to make clustering for tree saving.
	240	//if(! fEventClustered ){
	241	//fEvent->ExecuteClustering();
	242	//fEventClustered = true;
	243	//}
	244	fEvent->SetTime(fTime);
	245	fTTree->Fill();
	246	} else {
	247	if( fEvent ) delete fEvent;
	248	fEvent = event;
	249	fTTree->Fill();
	250	fEvent = 0;
	251	}
	252	}
	253	//-------------------------------------------------------------------
	254	void AliPHOSDApi0mip::FillHist(AliPHOSDATreeEvent* event){
	255	// Fill all information (AliPHOSDATreeEvent) into histograms
	256	// if event==NULL, the obtained information (=fEvent) is used.
	257
	258	if(! fCreateHist ) CreateHist();
	259	if( event == 0 ){
	260	if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
	261	if(! fEventClustered ){
	262	fEvent->ExecuteClustering();
	263	fEventClustered = true;
	264	}
	265	fEvent->SetTime(fTime);
	266	event = fEvent;
	267	}
	268
	269	// -- Constant
	270	int nclt1, nclt2;
	271	int ndigits;
	272	int ndigitsall = 0;
	273	AliPHOSDATreeDigit digit1, digit2;
	274	const float kDistanceFromVertex = 460.;
	275	float dist1, dist2, px, py;
	276	float mass, cosine, pt, asym;
	277
	278	// -- Filling into Time Histogram
	279	double start = fH1Time->GetBinContent(1);
	280	double stop = fH1Time->GetBinContent(2);
	281	double currenttime = (double) (event->GetTime());
	282	if( start == 0 && stop == 0 ) start = stop = currenttime;
	283	if( currenttime < start ) start = currenttime;
	284	if( currenttime > stop ) stop = currenttime;
	285	fH1Time->SetBinContent(1,start);
	286	fH1Time->SetBinContent(2,stop);
	287
288	// -- Start looping clusters
289	nclt1 = event->GetNClusters();
290	fH1ClusterNum->Fill(nclt1);
291	//std::cout<<" AliPHOSDApi0mip::FillHisto() DEBUG:: Number of clusters: "<<nclt1<<std::endl;
292	while( nclt1-- ){
293	AliPHOSDATreeCluster& clt1 = event->GetCluster(nclt1);
294	ndigits = clt1.GetNDigits();
295	while( ndigits-- ){
296	AliPHOSDATreeDigit& digit = clt1.GetDigit(ndigits);
297	fH2EneDigitId->Fill(digit.GetDigitId(),digit.GetEnergy());
298	ndigitsall++;
299	}
300	digit1 = clt1.GetMaxDigit();
301	if( clt1.GetEnergy() > 0 && clt1.GetNDigits() >= 2 ){
302	fH2MipDigitId->Fill(digit1.GetDigitId(),clt1.GetEnergy());
303	}
304	if( clt1.GetEnergy() > 0.2 ){
305	nclt2 = nclt1;
306	while( nclt2-- ){
307	AliPHOSDATreeCluster& clt2 = event->GetCluster(nclt2);
308	digit2 = clt2.GetMaxDigit();
309	if( clt2.GetEnergy() > 0.2 ){
310	//
311	// Following approximation is applied in order to avoid many calling "sqrt" function.
312	// sqrt( 1 + x^2 ) =~ 1 + 0.5x^2 (for x<<1)
313	//
314	//cosine = ((clt1.GetRow()-31.5)(clt2.GetRow()-31.5)4.84 + (clt1.GetCol()-27.5)(clt2.GetCol()-27.5)4.84 + kDistanceFromVertex*kDistanceFromVertex)
315	// / ( kDistanceFromVertex + (clt1.GetRow()-31.5)(clt1.GetRow()-31.5)2.42/kDistanceFromVertex + (clt1.GetCol()-27.5)(clt1.GetCol()-27.5)2.42/kDistanceFromVertex )
316	// / ( kDistanceFromVertex + (clt2.GetRow()-31.5)(clt2.GetRow()-31.5)2.42/kDistanceFromVertex + (clt2.GetCol()-27.5)(clt2.GetCol()-27.5)2.42/kDistanceFromVertex );
317	//
318
319
320	dist1 = sqrt( kDistanceFromVertexkDistanceFromVertex + (clt1.GetRow()-31.5)(clt1.GetRow()-31.5)4.84 + (clt1.GetCol()-27.5)(clt1.GetCol()-27.5)*4.84 );
321	dist2 = sqrt( kDistanceFromVertexkDistanceFromVertex + (clt2.GetRow()-31.5)(clt2.GetRow()-31.5)4.84 + (clt2.GetCol()-27.5)(clt2.GetCol()-27.5)*4.84 );
322	cosine = ((clt1.GetRow()-31.5)(clt2.GetRow()-31.5)4.84 + (clt1.GetCol()-27.5)(clt2.GetCol()-27.5)4.84 + kDistanceFromVertex*kDistanceFromVertex)
323	/ dist1 / dist2;
324	mass = sqrt( 2.0 * clt1.GetEnergy() * clt2.GetEnergy() * (1 - cosine) );
325	//std::cout<<" DEBUG: dist1:"<<dist1<<" dist2:"<<dist2<<" cosine:"<<cosine<<" c1:"<<clt1.GetEnergy()<<" c2:"<<clt2.GetEnergy()<<" mass:"<<mass<<std::endl;
326	//std::cout<<" : clt1(row,col)=("<<clt1.GetRow()<<","<<clt1.GetCol()<<") clt2(row,col)=("<<clt2.GetRow()<<","<<clt2.GetCol()<<") "<<std::endl;
327	asym = fabs(clt1.GetEnergy() - clt2.GetEnergy())/(clt1.GetEnergy() + clt2.GetEnergy());
328	py = clt1.GetEnergy() * kDistanceFromVertex / dist1 + clt2.GetEnergy() * kDistanceFromVertex / dist2;
329	px = clt1.GetEnergy() * (clt1.GetCol()-27.5)2.42/dist1 + clt2.GetEnergy() (clt2.GetCol()-27.5)*2.42/ dist2;
330	pt = sqrt(pxpx+pypy);
331	fH3Pi0AsymPt->Fill(asym,pt,mass);
332	if( digit1.IsValid() && digit1.IsValid() &&
333	digit1.GetDigitId() != digit2.GetDigitId() &&
334	//pt > 1.5 && clt1.GetEnergy() > 0.3 && clt2.GetEnergy() > 0.3 ){
335	//pt > 1.3 && clt1.GetEnergy() > 0.25 && clt2.GetEnergy() > 0.25 ){
336	pt > 1.0 && clt1.GetEnergy() > 0.2 && clt2.GetEnergy() > 0.2 ){
337	fH2Pi0DigitId->Fill(digit1.GetDigitId(),mass);
338	fH2Pi0DigitId->Fill(digit2.GetDigitId(),mass);
339	}
340	} // End of energy cut
341	} // End of second cluster loop (clt2)
342	} // End of energy cut
343	} // End of first cluster loop (clt1)
344	fH1DigitNum->Fill(ndigitsall);
345	//
346	}
347	//-------------------------------------------------------------------
c335b3c5	348	void AliPHOSDApi0mip::Print(Option_t *option) const
c335b3c5	349	{
12f6dd57	350	// Print Out
	351
	352	//fTFile->ls();
	353	//fTTree->Print();
	354	if( fEvent ){
	355	//fEvent->ExecuteClustering();
c335b3c5	356	fEvent->Print(option);
12f6dd57	357	}
	358	}
	359	//-------------------------------------------------------------------