[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGenerator.cxx

#include <iostream>

#include <TDatabasePDG.h>
#include <TRandom.h>
#include <TH2F.h>
#include <TGeoGlobalMagField.h>

#include <AliLog.h>
#include <AliTPCROC.h>
#include <AliTrackPointArray.h>
#include <AliTrackerBase.h>
#include <AliCDBManager.h>
#include <AliTPCParam.h>
#include <AliGeomManager.h>
#include <AliTPCcalibDB.h>
#include <AliTPCclusterMI.h>
#include <AliTPCSpaceCharge3D.h>

#include "AliToyMCEvent.h"
#include "AliToyMCTrack.h"

#include "AliToyMCEventGenerator.h"

ClassImp(AliToyMCEventGenerator);


AliToyMCEventGenerator::AliToyMCEventGenerator()
  :TObject()
  ,fTPCParam(0x0)
  ,fEvent(0x0)
  ,fSpaceCharge(0x0)
  ,fOutputFileName("toyMC.root")
  ,fOutFile(0x0)
  ,fOutTree(0x0)
{
  //TODO: Add method to set custom space charge file
  fSpaceCharge = new AliTPCSpaceCharge3D();
  fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz.root");
  fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
  //fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
  fSpaceCharge->InitSpaceCharge3DDistortion();
//   fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
//   fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
  //!!! This should be handled by the CongiOCDB macro
  // const char* ocdb="local://$ALICE_ROOT/OCDB/";
  // AliCDBManager::Instance()->SetDefaultStorage(ocdb);
  // AliCDBManager::Instance()->SetRun(0);   
  // TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  // AliGeomManager::LoadGeometry("");
  fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  //std::cout<<"----------->"<<fTPCParam<<std::endl;
  fTPCParam->ReadGeoMatrices();
  //std::cout<<"----------->"<<fTPCParam<<std::endl;

}
//________________________________________________________________
AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
  :TObject(gen)
  ,fTPCParam(gen.fTPCParam)
  ,fEvent(0x0)
  ,fSpaceCharge(gen.fSpaceCharge)
  ,fOutputFileName(gen.fOutputFileName)
  ,fOutFile(0x0)
  ,fOutTree(0x0)
{
  //
}
//________________________________________________________________
AliToyMCEventGenerator::~AliToyMCEventGenerator() 
{
  delete fSpaceCharge;
  //!!! Is fTPCParam not owned by the CDB manager?
  delete fTPCParam;
}
//________________________________________________________________
Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
{
  //
  //
  //

  //!!! Should this complete function not be moved to AliToyMCTrack, and directly called from the constructor
  //    or setter of the track parameters?

  if(!fTPCParam) {
    fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
    fTPCParam->ReadGeoMatrices();
    std::cout << "init tpc param" << std::endl;
  }
  const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
  // const Double_t kRTPC0  =AliTPCROC::Instance()->GetPadRowRadii(0,0);
  // const Double_t kRTPC1  =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
  const Double_t kMaxSnp = 0.85;  
  const Double_t kSigmaY=0.1; 
  const Double_t kSigmaZ=0.1;
  //!!! why divide by 1000000? This will be in cm/us then. But I guess everything should be in sec;
  //   const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
  const Double_t kDriftVel = fTPCParam->GetDriftV();
  const Double_t kMaxZ0=fTPCParam->GetZLength();
  const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
    
  const Double_t iFCRadius=  83.5; //radius constants found in AliTPCCorrection.cxx 
  const Double_t oFCRadius= 254.5;
  
  AliToyMCTrack track(trackIn);
  const Int_t nMaxPoints = kNRows+50;
  AliTrackPointArray pointArray0(nMaxPoints);
  AliTrackPointArray pointArray1(nMaxPoints);
  Double_t xyz[3];

  //!!! when does the propagation not work, how often does it happen?
  if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) {
    AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
    return kFALSE;
  }

  Int_t npoints=0;
  Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame


  //Simulate track from inner field cage radius to outer and save points along trajectory
  Int_t nMinPoints = 40;
  //!!! changed to loop over pad rows
  //      an increment of 1cm only might loose pad rows in IROC: 4x7.5 mm2
  //      or we change this to a smaller increment
  //    Do we need this initial creation of space points at all? can't we
  //      simply loop over all rows, crate the clusters in the center of the row,
  //      distort them and fill the distorted clusters
  
  //for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){
  for (Int_t irow=0; irow<159; ++irow ){
    const Int_t isec    = irow<63?0:36;
    const Int_t isecrow = irow<63?irow:irow-63;
    Double_t radius = fTPCParam->GetPadRowRadii(isec,isecrow);

    //!!! changed from return 0 to continue -> Please check
    if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) {
      AliError(Form("Propagation to %d: %.2f failed\n",irow,radius));
      continue;
    }
    track.GetXYZ(xyz);

    //!!! Why is this smeared
    xyz[0]+=gRandom->Gaus(0,0.000005);
    xyz[1]+=gRandom->Gaus(0,0.000005);
    xyz[2]+=gRandom->Gaus(0,0.000005);

    if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
    //!!! for future we might want to incude space points outside the active area
    //      since due to the distotions they might get inside
    //      in this case it will of course make sense to make a dense
    //      space point creation first and afterwards calculate the cluster
    //      but then we need to treat cluster and distorted cluster symmetrically (see below)
    if (TMath::Abs(track.GetX())<iFCRadius) continue;
    if (TMath::Abs(track.GetX())>oFCRadius) continue;

    AliTrackPoint pIn0;                               // space point          
    AliTrackPoint pIn1;
    Int_t sector= (xyz[2]>0)? 0:18;
    pointArray0.GetPoint(pIn0,npoints);
    pointArray1.GetPoint(pIn1,npoints);
    Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
    Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
    fSpaceCharge->DistortPoint(distPoint, sector);
    pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
    pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
    track.Rotate(alpha);
    AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
    AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
    prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
    prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
    pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
    pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
    pointArray0.AddPoint(npoints, &pIn0);      
    pointArray1.AddPoint(npoints, &pIn1);
    npoints++;
    if (npoints>=nMaxPoints) break;
  }

  if (npoints<nMinPoints) return 0;

  //save space points and make clusters of distorted points

  //!!! Here both, clusters and distorted clusters should be in the center of the pad row
  //      So in case the distorted clusters are calculated as the COG of all space points
  //      in a row, the same should be done for the undistorted clusters! There should
  //      only be one cluster per track per row
  Int_t lastRow = 0;
  Int_t pntsInCurrentRow = 0;
  Double_t xt=0, yt=0, zt=0;
  for(Int_t iPoint = 0; iPoint<npoints; iPoint++){

    //undistorted cluster, currently we assume one space point per pad row, see comments above
    AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());

    const Float_t *xArr = pointArray0.GetX();
    const Float_t *yArr = pointArray0.GetY();
    const Float_t *zArr = pointArray0.GetZ();

    tempCl->SetX(xArr[iPoint]);
    tempCl->SetY(yArr[iPoint]);
    tempCl->SetZ(zArr[iPoint]);
    
    Float_t xyzUC[3] = {xArr[iPoint],yArr[iPoint],zArr[iPoint]};
    Int_t indexUC[3] = {0};
    
    Float_t padRowUC = fTPCParam->GetPadRow(xyzUC,indexUC);
    Int_t   sectorUC = indexUC[1];

    // use pad row = 255 if outside active area
    if(padRowUC < 0 || (sectorUC < 36 && padRowUC >62) || (sectorUC > 35 &&padRowUC > 95) ) padRowUC=255; //outside sensitive area

    tempCl->SetPad(xyzUC[1]);
    tempCl->SetRow(padRowUC);
    tempCl->SetDetector(sectorUC);
    tempCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(tempCl->GetZ()))/kDriftVel); // set time as t0  + drift time from dist z

    // Distorted cluster
    const Float_t *xArrDist = pointArray1.GetX();
    const Float_t *yArrDist = pointArray1.GetY();
    const Float_t *zArrDist = pointArray1.GetZ();
    
    Float_t xyz2[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
    Int_t index[3] = {0};

    Float_t padRow = fTPCParam->GetPadRow(xyz2,index);
    Int_t sector = index[1];

    if(padRow < 0 || (sector < 36 && padRow>62) || (sector > 35 &&padRow > 95) ) continue; //outside sensitive area

    //!!! could it be that the last row is missed?
    //    assume row=95 and there is only 1 space point then this part will not be reached I think
    if(lastRow!=padRow) {
      

      //make and save distorted cluster
      if(pntsInCurrentRow>0){
	xt/=pntsInCurrentRow;
	yt/=pntsInCurrentRow;
	zt/=pntsInCurrentRow;
	AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
	tempDistCl->SetX(xt);
	tempDistCl->SetY(yt);
	tempDistCl->SetZ(zt);
	Float_t clxyz[3] = {xt,yt,zt};
	Int_t clindex[3] = {0};
	Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
// 	Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
// 	Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
        //!!! pad should be fractional
// 	tempDistCl->SetPad(pad);
        tempDistCl->SetPad(clxyz[1]);
        tempDistCl->SetRow(clRow);
	tempDistCl->SetDetector(clindex[1]);
	tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0  + drift time from dist z coordinate to pad plane
      }
      xt = 0;
      yt = 0;
      zt = 0;
      pntsInCurrentRow = 0;
    }

    xt+=xArrDist[iPoint];
    yt+=yArrDist[iPoint];
    zt+=zArrDist[iPoint];
    pntsInCurrentRow++;
    lastRow = padRow;
  }
 
  return 1;

  
}
//________________________________________________________________
Bool_t AliToyMCEventGenerator::ConnectOutputFile()
{
  //
  // Create the output file name and tree and connect the event
  //

  fOutFile = new TFile(fOutputFileName.Data(),"recreate");

  if (!fOutFile || !fOutFile->IsOpen()){
    delete fOutFile;
    fOutFile=0x0;
    return kFALSE;
  }
  
  fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
  fOutTree->Branch("event","AliToyMCEvent",&fEvent);

  return kTRUE;
}

//________________________________________________________________
Bool_t AliToyMCEventGenerator::CloseOutputFile()
{
  //
  // close the output file
  //
  if (!fOutFile) return kFALSE;
  fOutFile->Write();
  fOutFile->Close();
  delete fOutFile;
  fOutFile=0x0;

  return kTRUE;
}

//________________________________________________________________
void AliToyMCEventGenerator::FillTree()
{
  // fill the tree
  if (fOutTree) fOutTree->Fill();
}
Commit	Line	Data
526ddf0e	1	#include <iostream>
a1a695e5	2
526ddf0e	3	#include <TDatabasePDG.h>
	4	#include <TRandom.h>
	5	#include <TH2F.h>
a1a695e5	6	#include <TGeoGlobalMagField.h>
a1a695e5	7
b0c26b55	8	#include <AliLog.h>
a1a695e5	9	#include <AliTPCROC.h>
a1a695e5	10	#include <AliTrackPointArray.h>
526ddf0e	11	#include <AliTrackerBase.h>
	12	#include <AliCDBManager.h>
	13	#include <AliTPCParam.h>
	14	#include <AliGeomManager.h>
	15	#include <AliTPCcalibDB.h>
526ddf0e	16	#include <AliTPCclusterMI.h>
a1a695e5	17	#include <AliTPCSpaceCharge3D.h>
	18
	19	#include "AliToyMCEvent.h"
	20	#include "AliToyMCTrack.h"
	21
	22	#include "AliToyMCEventGenerator.h"
	23
de0014b7	24	ClassImp(AliToyMCEventGenerator);
526ddf0e	25
526ddf0e	26
de0014b7	27	AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e	28	:TObject()
526ddf0e	29	,fTPCParam(0x0)
a1a695e5	30	,fEvent(0x0)
	31	,fSpaceCharge(0x0)
	32	,fOutputFileName("toyMC.root")
	33	,fOutFile(0x0)
	34	,fOutTree(0x0)
526ddf0e	35	{
	36	//TODO: Add method to set custom space charge file
	37	fSpaceCharge = new AliTPCSpaceCharge3D();
a1a695e5	38	fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz.root");
526ddf0e	39	fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
	40	//fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
	41	fSpaceCharge->InitSpaceCharge3DDistortion();
b0c26b55	42	// fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
b0c26b55	43	// fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
a1a695e5	44	//!!! This should be handled by the CongiOCDB macro
efef400b	45	// const char* ocdb="local://$ALICE_ROOT/OCDB/";
	46	// AliCDBManager::Instance()->SetDefaultStorage(ocdb);
	47	// AliCDBManager::Instance()->SetRun(0);
	48	// TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
	49	// AliGeomManager::LoadGeometry("");
526ddf0e	50	fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
efef400b	51	//std::cout<<"----------->"<<fTPCParam<<std::endl;
526ddf0e	52	fTPCParam->ReadGeoMatrices();
efef400b	53	//std::cout<<"----------->"<<fTPCParam<<std::endl;
efef400b	54
526ddf0e	55	}
526ddf0e	56	//________________________________________________________________
de0014b7	57	AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e	58	:TObject(gen)
a1a695e5	59	,fTPCParam(gen.fTPCParam)
a1a695e5	60	,fEvent(0x0)
526ddf0e	61	,fSpaceCharge(gen.fSpaceCharge)
a1a695e5	62	,fOutputFileName(gen.fOutputFileName)
	63	,fOutFile(0x0)
	64	,fOutTree(0x0)
526ddf0e	65	{
	66	//
	67	}
	68	//________________________________________________________________
de0014b7	69	AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e	70	{
526ddf0e	71	delete fSpaceCharge;
a1a695e5	72	//!!! Is fTPCParam not owned by the CDB manager?
526ddf0e	73	delete fTPCParam;
	74	}
	75	//________________________________________________________________
a1a695e5	76	Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
	77	{
	78	//
	79	//
	80	//
	81
	82	//!!! Should this complete function not be moved to AliToyMCTrack, and directly called from the constructor
	83	// or setter of the track parameters?
526ddf0e	84
	85	if(!fTPCParam) {
	86	fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
	87	fTPCParam->ReadGeoMatrices();
	88	std::cout << "init tpc param" << std::endl;
	89	}
	90	const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
	91	// const Double_t kRTPC0 =AliTPCROC::Instance()->GetPadRowRadii(0,0);
	92	// const Double_t kRTPC1 =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
	93	const Double_t kMaxSnp = 0.85;
b0c26b55	94	const Double_t kSigmaY=0.1;
526ddf0e	95	const Double_t kSigmaZ=0.1;
b0c26b55	96	//!!! why divide by 1000000? This will be in cm/us then. But I guess everything should be in sec;
	97	// const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
	98	const Double_t kDriftVel = fTPCParam->GetDriftV();
526ddf0e	99	const Double_t kMaxZ0=fTPCParam->GetZLength();
	100	const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
	101
	102	const Double_t iFCRadius= 83.5; //radius constants found in AliTPCCorrection.cxx
	103	const Double_t oFCRadius= 254.5;
	104
de0014b7	105	AliToyMCTrack track(trackIn);
526ddf0e	106	const Int_t nMaxPoints = kNRows+50;
	107	AliTrackPointArray pointArray0(nMaxPoints);
	108	AliTrackPointArray pointArray1(nMaxPoints);
	109	Double_t xyz[3];
	110
a1a695e5	111	//!!! when does the propagation not work, how often does it happen?
b0c26b55	112	if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) {
	113	AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
	114	return kFALSE;
	115	}
526ddf0e	116
	117	Int_t npoints=0;
	118	Float_t covPoint[6]={0,0,0, kSigmaYkSigmaY,0,kSigmaZkSigmaZ}; //covariance at the local frame
	119
	120
	121	//Simulate track from inner field cage radius to outer and save points along trajectory
b0c26b55	122	Int_t nMinPoints = 40;
	123	//!!! changed to loop over pad rows
	124	// an increment of 1cm only might loose pad rows in IROC: 4x7.5 mm2
	125	// or we change this to a smaller increment
	126	// Do we need this initial creation of space points at all? can't we
	127	// simply loop over all rows, crate the clusters in the center of the row,
	128	// distort them and fill the distorted clusters
	129
	130	//for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){
	131	for (Int_t irow=0; irow<159; ++irow ){
	132	const Int_t isec = irow<63?0:36;
	133	const Int_t isecrow = irow<63?irow:irow-63;
	134	Double_t radius = fTPCParam->GetPadRowRadii(isec,isecrow);
526ddf0e	135
b0c26b55	136	//!!! changed from return 0 to continue -> Please check
	137	if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) {
	138	AliError(Form("Propagation to %d: %.2f failed\n",irow,radius));
	139	continue;
	140	}
526ddf0e	141	track.GetXYZ(xyz);
b0c26b55	142
b0c26b55	143	//!!! Why is this smeared
526ddf0e	144	xyz[0]+=gRandom->Gaus(0,0.000005);
	145	xyz[1]+=gRandom->Gaus(0,0.000005);
	146	xyz[2]+=gRandom->Gaus(0,0.000005);
	147
	148	if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
b0c26b55	149	//!!! for future we might want to incude space points outside the active area
	150	// since due to the distotions they might get inside
	151	// in this case it will of course make sense to make a dense
	152	// space point creation first and afterwards calculate the cluster
	153	// but then we need to treat cluster and distorted cluster symmetrically (see below)
526ddf0e	154	if (TMath::Abs(track.GetX())<iFCRadius) continue;
	155	if (TMath::Abs(track.GetX())>oFCRadius) continue;
	156
	157	AliTrackPoint pIn0; // space point
	158	AliTrackPoint pIn1;
	159	Int_t sector= (xyz[2]>0)? 0:18;
	160	pointArray0.GetPoint(pIn0,npoints);
	161	pointArray1.GetPoint(pIn1,npoints);
	162	Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
	163	Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
	164	fSpaceCharge->DistortPoint(distPoint, sector);
	165	pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
	166	pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
	167	track.Rotate(alpha);
	168	AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
	169	AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
	170	prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
	171	prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
	172	pIn0=prot0.Rotate(-alpha); // rotate back to global frame
	173	pIn1=prot1.Rotate(-alpha); // rotate back to global frame
	174	pointArray0.AddPoint(npoints, &pIn0);
	175	pointArray1.AddPoint(npoints, &pIn1);
	176	npoints++;
	177	if (npoints>=nMaxPoints) break;
	178	}
	179
	180	if (npoints<nMinPoints) return 0;
	181
	182	//save space points and make clusters of distorted points
b0c26b55	183
	184	//!!! Here both, clusters and distorted clusters should be in the center of the pad row
	185	// So in case the distorted clusters are calculated as the COG of all space points
	186	// in a row, the same should be done for the undistorted clusters! There should
	187	// only be one cluster per track per row
526ddf0e	188	Int_t lastRow = 0;
	189	Int_t pntsInCurrentRow = 0;
	190	Double_t xt=0, yt=0, zt=0;
	191	for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
b0c26b55	192
b0c26b55	193	//undistorted cluster, currently we assume one space point per pad row, see comments above
526ddf0e	194	AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());
	195
	196	const Float_t *xArr = pointArray0.GetX();
	197	const Float_t *yArr = pointArray0.GetY();
	198	const Float_t *zArr = pointArray0.GetZ();
	199
526ddf0e	200	tempCl->SetX(xArr[iPoint]);
	201	tempCl->SetY(yArr[iPoint]);
	202	tempCl->SetZ(zArr[iPoint]);
	203
b0c26b55	204	Float_t xyzUC[3] = {xArr[iPoint],yArr[iPoint],zArr[iPoint]};
	205	Int_t indexUC[3] = {0};
	206
	207	Float_t padRowUC = fTPCParam->GetPadRow(xyzUC,indexUC);
	208	Int_t sectorUC = indexUC[1];
	209
	210	// use pad row = 255 if outside active area
	211	if(padRowUC < 0 \|\| (sectorUC < 36 && padRowUC >62) \|\| (sectorUC > 35 &&padRowUC > 95) ) padRowUC=255; //outside sensitive area
	212
	213	tempCl->SetPad(xyzUC[1]);
	214	tempCl->SetRow(padRowUC);
	215	tempCl->SetDetector(sectorUC);
	216	tempCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(tempCl->GetZ()))/kDriftVel); // set time as t0 + drift time from dist z
	217
	218	// Distorted cluster
	219	const Float_t *xArrDist = pointArray1.GetX();
	220	const Float_t *yArrDist = pointArray1.GetY();
	221	const Float_t *zArrDist = pointArray1.GetZ();
	222
a1a695e5	223	Float_t xyz2[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
526ddf0e	224	Int_t index[3] = {0};
526ddf0e	225
a1a695e5	226	Float_t padRow = fTPCParam->GetPadRow(xyz2,index);
526ddf0e	227	Int_t sector = index[1];
	228
	229	if(padRow < 0 \|\| (sector < 36 && padRow>62) \|\| (sector > 35 &&padRow > 95) ) continue; //outside sensitive area
	230
b0c26b55	231	//!!! could it be that the last row is missed?
b0c26b55	232	// assume row=95 and there is only 1 space point then this part will not be reached I think
526ddf0e	233	if(lastRow!=padRow) {
	234
	235
	236	//make and save distorted cluster
	237	if(pntsInCurrentRow>0){
	238	xt/=pntsInCurrentRow;
	239	yt/=pntsInCurrentRow;
	240	zt/=pntsInCurrentRow;
	241	AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
	242	tempDistCl->SetX(xt);
	243	tempDistCl->SetY(yt);
	244	tempDistCl->SetZ(zt);
	245	Float_t clxyz[3] = {xt,yt,zt};
	246	Int_t clindex[3] = {0};
	247	Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
b0c26b55	248	// Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
	249	// Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
	250	//!!! pad should be fractional
	251	// tempDistCl->SetPad(pad);
	252	tempDistCl->SetPad(clxyz[1]);
	253	tempDistCl->SetRow(clRow);
526ddf0e	254	tempDistCl->SetDetector(clindex[1]);
	255	tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0 + drift time from dist z coordinate to pad plane
	256	}
	257	xt = 0;
	258	yt = 0;
	259	zt = 0;
	260	pntsInCurrentRow = 0;
	261	}
	262
	263	xt+=xArrDist[iPoint];
	264	yt+=yArrDist[iPoint];
	265	zt+=zArrDist[iPoint];
	266	pntsInCurrentRow++;
	267	lastRow = padRow;
	268	}
	269
	270	return 1;
	271
	272
	273
	274	}
a1a695e5	275	//________________________________________________________________
	276	Bool_t AliToyMCEventGenerator::ConnectOutputFile()
	277	{
	278	//
	279	// Create the output file name and tree and connect the event
	280	//
	281
	282	fOutFile = new TFile(fOutputFileName.Data(),"recreate");
	283
	284	if (!fOutFile \|\| !fOutFile->IsOpen()){
	285	delete fOutFile;
	286	fOutFile=0x0;
	287	return kFALSE;
	288	}
	289
	290	fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
	291	fOutTree->Branch("event","AliToyMCEvent",&fEvent);
	292
	293	return kTRUE;
	294	}
	295
	296	//________________________________________________________________
	297	Bool_t AliToyMCEventGenerator::CloseOutputFile()
	298	{
	299	//
	300	// close the output file
	301	//
	302	if (!fOutFile) return kFALSE;
	303	fOutFile->Write();
	304	fOutFile->Close();
	305	delete fOutFile;
	306	fOutFile=0x0;
	307
	308	return kTRUE;
	309	}
	310
	311	//________________________________________________________________
	312	void AliToyMCEventGenerator::FillTree()
	313	{
	314	// fill the tree
	315	if (fOutTree) fOutTree->Fill();
	316	}
	317