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

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  fast TPC cluster simulation                                              //
//                                                                          //
///////////////////////////////////////////////////////////////////////////////

#include <TParticle.h>
#include <TVector.h>

#include "AliRunLoader.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliTPC.h"
#include "AliTPCParam.h"
#include "AliTPCClustersArray.h"
#include "AliTPCClustersRow.h"
#include "AliTPCcluster.h"
#include "AliComplexCluster.h"
#include "AliTPCFast.h"
#include "AliLog.h"

ClassImp(AliTPCFast)


//_____________________________________________________________________________
void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
{
  //--------------------------------------------------------
  // TPC simple cluster generator from hits
  // obtained from the TPC Fast Simulator
  // The point errors are taken from the parametrization
  //--------------------------------------------------------

  //-----------------------------------------------------------------
  // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
  //-----------------------------------------------------------------
  // Adopted to Marian's cluster data structure by I.Belikov, CERN,
  // Jouri.Belikov@cern.ch
  //----------------------------------------------------------------
  
  /////////////////////////////////////////////////////////////////////////////
  //
  //---------------------------------------------------------------------
  //   ALICE TPC Cluster Parameters
  //--------------------------------------------------------------------
       
  
  // Cluster width in rphi
  const Float_t kACrphi=0.18322;
  const Float_t kBCrphi=0.59551e-3;
  const Float_t kCCrphi=0.60952e-1;
  // Cluster width in z
  const Float_t kACz=0.19081;
  const Float_t kBCz=0.55938e-3;
  const Float_t kCCz=0.30428;


  AliLoader* loader = runLoader->GetLoader("TPCLoader");
  if (!loader) {
    AliError("No TPC loader found");
    return;
  }
  if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  AliRun* aliRun = runLoader->GetAliRun();
  if (!aliRun) {
    AliError("Couldn't get AliRun object");
    return;
  }
  AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
  if (!tpc) {
    AliError("Couldn't get TPC detector");
    return;
  }
  AliTPCParam* param = tpc->GetParam();
  if (!param) {
    AliError("No TPC parameters available");
    return;
  }

  //if(fDefaults == 0) SetDefaults();

  Float_t sigmaRphi,sigmaZ,clRphi,clZ;
  //
  TParticle *particle; // pointer to a given particle
  AliTPChit *tpcHit; // pointer to a sigle TPC hit
  Int_t sector;
  Int_t ipart;
  Float_t xyz[5];
  Float_t pl,pt,tanth,rpad,ratio;
  Float_t cph,sph;
  
  //---------------------------------------------------------------
  //  Get the access to the tracks 
  //---------------------------------------------------------------
  
  TTree *tH = loader->TreeH();
  if (tH == 0x0)
    AliFatal("Can not find TreeH in folder");

  tpc->SetTreeAddress();
  
  Stat_t ntracks = tH->GetEntries();

  //Switch to the output file
  
  if (loader->TreeR() == 0x0) loader->MakeTree("R");
  
  AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
  
  runLoader->CdGAFile();
  //param->Write(param->GetTitle());

  AliTPCClustersArray carray;
  carray.Setup(param);
  carray.SetClusterType("AliTPCcluster");
  carray.MakeTree(loader->TreeR());

  Int_t nclusters=0; //cluster counter
  
  //------------------------------------------------------------
  // Loop over all sectors (72 sectors for 20 deg
  // segmentation for both lower and upper sectors)
  // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
  // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
  //
  // First cluster for sector 0 starts at "0"
  //------------------------------------------------------------
   
  for(Int_t isec=0;isec<param->GetNSector();isec++){
    //MI change
    param->AdjustCosSin(isec,cph,sph);
    
    //------------------------------------------------------------
    // Loop over tracks
    //------------------------------------------------------------
    
    for(Int_t track=0;track<ntracks;track++){
      tpc->ResetHits();
      tH->GetEvent(track);
      //
      //  Get number of the TPC hits
      //     
       tpcHit = (AliTPChit*)tpc->FirstHit(-1);

      // Loop over hits
      //
       while(tpcHit){
 
	 if (tpcHit->fQ == 0.) {
           tpcHit = (AliTPChit*) tpc->NextHit();
           continue; //information about track (I.Belikov)
	 }
	sector=tpcHit->fSector; // sector number

       if(sector != isec){
	 tpcHit = (AliTPChit*) tpc->NextHit();
	 continue; 
       }
	ipart=tpcHit->Track();
	particle=aliRun->GetMCApp()->Particle(ipart);
	pl=particle->Pz();
	pt=particle->Pt();
	if(pt < 1.e-9) pt=1.e-9;
	tanth=pl/pt;
	tanth = TMath::Abs(tanth);
	rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
	ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason

	//   space-point resolutions
	
	sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
	sigmaZ   =AliTPCcluster::SigmaZ2(rpad,tanth   );
	
	//   cluster widths
	
	clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
	clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
	
	// temporary protection
	
	if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
	if(sigmaZ < 0.) sigmaZ=0.4e-3;
	if(clRphi < 0.) clRphi=2.5e-3;
	if(clZ < 0.) clZ=2.5e-5;
	
	//
	
	//
	// smearing --> rotate to the 1 (13) or to the 25 (49) sector,
	// then the inaccuracy in a X-Y plane is only along Y (pad row)!
	//
        Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
	Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
	xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
          Float_t alpha=(isec < param->GetNInnerSector()) ?
	  param->GetInnerAngle() : param->GetOuterAngle();
          Float_t ymax=xprim*TMath::Tan(0.5*alpha);
          if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim; 
	xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
          if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z(); 
	xyz[2]=sigmaRphi;                                     // fSigmaY2
	xyz[3]=sigmaZ;                                        // fSigmaZ2
	xyz[4]=tpcHit->fQ;                                    // q

        AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
        if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);	

        Int_t tracks[3]={tpcHit->Track(), -1, -1};
	AliTPCcluster cluster(tracks,xyz);

        clrow->InsertCluster(&cluster); nclusters++;

        tpcHit = (AliTPChit*)tpc->NextHit();
        

      } // end of loop over hits

    }   // end of loop over tracks

    Int_t nrows=param->GetNRow(isec);
    for (Int_t irow=0; irow<nrows; irow++) {
        AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
        if (!clrow) continue;
        carray.StoreRow(isec,irow);
        carray.ClearRow(isec,irow);
    }

  } // end of loop over sectors  

  //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
  loader->WriteRecPoints("OVERWRITE");
  
  
} // end of function

//_________________________________________________________________
void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
					  AliTPCClustersArray* clustersArray,
					  Int_t isec) const
{
  //--------------------------------------------------------
  //calculate exact cross point of track and given pad row
  //resulting values are expressed in "digit" coordinata
  //--------------------------------------------------------

  //-----------------------------------------------------------------
  // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
  //-----------------------------------------------------------------
  //
  if (clustersArray==0){    
    return;
  }
  AliLoader* loader = runLoader->GetLoader("TPCLoader");
  if (!loader) {
    AliError("No TPC loader found");
    return;
  }
  if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  AliRun* aliRun = runLoader->GetAliRun();
  if (!aliRun) {
    AliError("Couldn't get AliRun object");
    return;
  }
  AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
  if (!tpc) {
    AliError("Couldn't get TPC detector");
    return;
  }
  AliTPCParam* param = tpc->GetParam();
  if (!param) {
    AliError("No TPC parameters available");
    return;
  }
  //
  TParticle *particle; // pointer to a given particle
  AliTPChit *tpcHit; // pointer to a sigle TPC hit
  //  Int_t sector,nhits;
  Int_t ipart;
  const Int_t kcmaxhits=30000;
  TVector * xxxx = new TVector(kcmaxhits*4);
  TVector & xxx = *xxxx;
  Int_t maxhits = kcmaxhits;
  //construct array for each padrow
  for (Int_t i=0; i<param->GetNRow(isec);i++) 
    clustersArray->CreateRow(isec,i);
  
  //---------------------------------------------------------------
  //  Get the access to the tracks 
  //---------------------------------------------------------------
  
  TTree *tH = loader->TreeH();
  if (tH == 0x0)
    AliFatal("Can not find TreeH in folder");

  tpc->SetTreeAddress();

  Stat_t ntracks = tH->GetEntries();
  Int_t npart = aliRun->GetMCApp()->GetNtrack();
  //MI change
  TBranch * branch=0;
  if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
  else branch = tH->GetBranch("TPC");

  //------------------------------------------------------------
  // Loop over tracks
  //------------------------------------------------------------

  for(Int_t track=0;track<ntracks;track++){ 
    Bool_t isInSector=kTRUE;
    tpc->ResetHits();
     isInSector = tpc->TrackInVolume(isec,track);
    if (!isInSector) continue;
    //MI change
    branch->GetEntry(track); // get next track
    //
    //  Get number of the TPC hits and a pointer
    //  to the particles
    // Loop over hits
    //
    Int_t currentIndex=0;
    Int_t lastrow=-1;  //last writen row

    //M.I. changes

    tpcHit = (AliTPChit*)tpc->FirstHit(-1);
    while(tpcHit){
      
      Int_t sector=tpcHit->fSector; // sector number
      if(sector != isec){
	tpcHit = (AliTPChit*) tpc->NextHit();
	continue; 
      }

      ipart=tpcHit->Track();
      if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
      
      //find row number

      Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
      Int_t    index[3]={1,isec,0};
      Int_t    currentrow = param->GetPadRow(x,index) ;	
      if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
      if (lastrow<0) lastrow=currentrow;
      if (currentrow==lastrow){
	if ( currentIndex>=maxhits){
	  maxhits+=kcmaxhits;
	  xxx.ResizeTo(4*maxhits);
	}     
	xxx(currentIndex*4)=x[0];
	xxx(currentIndex*4+1)=x[1];
	xxx(currentIndex*4+2)=x[2];	
	xxx(currentIndex*4+3)=tpcHit->fQ;
	currentIndex++;	
      }
      else 
	if (currentIndex>2){
	  Float_t sumx=0;
	  Float_t sumx2=0;
	  Float_t sumx3=0;
	  Float_t sumx4=0;
	  Float_t sumy=0;
	  Float_t sumxy=0;
	  Float_t sumx2y=0;
	  Float_t sumz=0;
	  Float_t sumxz=0;
	  Float_t sumx2z=0;
	  Float_t sumq=0;
	  for (Int_t index=0;index<currentIndex;index++){
	    Float_t x,x2,x3,x4;
	    x=x2=x3=x4=xxx(index*4);
	    x2*=x;
	    x3*=x2;
	    x4*=x3;
	    sumx+=x;
	    sumx2+=x2;
	    sumx3+=x3;
	    sumx4+=x4;
	    sumy+=xxx(index*4+1);
	    sumxy+=xxx(index*4+1)*x;
	    sumx2y+=xxx(index*4+1)*x2;
	    sumz+=xxx(index*4+2);
	    sumxz+=xxx(index*4+2)*x;
	    sumx2z+=xxx(index*4+2)*x2;	 
	    sumq+=xxx(index*4+3);
	  }
	  Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
	  Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
	    sumx2*(sumx*sumx3-sumx2*sumx2);
	  
	  Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
	    sumx2*(sumxy*sumx3-sumx2y*sumx2);
	  Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
	    sumx2*(sumxz*sumx3-sumx2z*sumx2);
	  
	  Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
	    sumx2*(sumx*sumx2y-sumx2*sumxy);
	  Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
	    sumx2*(sumx*sumx2z-sumx2*sumxz);
	  
	  if (TMath::Abs(det)<0.00001){
	     tpcHit = (AliTPChit*)tpc->NextHit();
	    continue;
	  }
	
	  Float_t y=detay/det+centralPad;
	  Float_t z=detaz/det;	
	  Float_t by=detby/det; //y angle
	  Float_t bz=detbz/det; //z angle
	  sumy/=Float_t(currentIndex);
	  sumz/=Float_t(currentIndex);

	  AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
	  if (row!=0) {
	    AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
	    //	  AliComplexCluster cl;
	    cl->fX=z;
	    cl->fY=y;
	    cl->fQ=sumq;
	    cl->fSigmaX2=bz;
	    cl->fSigmaY2=by;
	    cl->fTracks[0]=ipart;
	  }
	  currentIndex=0;
	  lastrow=currentrow;
	} //end of calculating cluster for given row
	
	
      tpcHit = (AliTPChit*)tpc->NextHit();
    } // end of loop over hits
  }   // end of loop over tracks 
  //write padrows to tree 
  for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
    clustersArray->StoreRow(isec,ii);    
    clustersArray->ClearRow(isec,ii);        
  }
  xxxx->Delete();
 
}


void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
					   AliTPCClustersArray* clustersArray) const
{

  //
  //fill clones array with intersection of current point with the
  //middle of the row
  AliLoader* loader = runLoader->GetLoader("TPCLoader");
  if (!loader) {
    AliError("No TPC loader found");
    return;
  }
  if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  AliRun* aliRun = runLoader->GetAliRun();
  if (!aliRun) {
    AliError("Couldn't get AliRun object");
    return;
  }
  AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
  if (!tpc) {
    AliError("Couldn't get TPC detector");
    return;
  }
  AliTPCParam* param = tpc->GetParam();
  if (!param) {
    AliError("No TPC parameters available");
    return;
  }

  Int_t sector;
  Int_t ipart;
  
  const Int_t kcmaxhits=30000;
  TVector * xxxx = new TVector(kcmaxhits*4);
  TVector & xxx = *xxxx;
  Int_t maxhits = kcmaxhits;
      
  //
  AliTPChit * tpcHit=0;
  tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
  Int_t currentIndex=0;
  Int_t lastrow=-1;  //last writen row

  while (tpcHit){
    if (tpcHit==0) continue;
    sector=tpcHit->fSector; // sector number
    ipart=tpcHit->Track();
    
    //find row number
    
    Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
    Int_t    index[3]={1,sector,0};
    Int_t    currentrow = param->GetPadRow(x,index) ;	
    if (currentrow<0) continue;
    if (lastrow<0) lastrow=currentrow;
    if (currentrow==lastrow){
      if ( currentIndex>=maxhits){
	maxhits+=kcmaxhits;
	xxx.ResizeTo(4*maxhits);
      }     
      xxx(currentIndex*4)=x[0];
      xxx(currentIndex*4+1)=x[1];
      xxx(currentIndex*4+2)=x[2];	
      xxx(currentIndex*4+3)=tpcHit->fQ;
      currentIndex++;	
    }
    else 
      if (currentIndex>2){
	Float_t sumx=0;
	Float_t sumx2=0;
	Float_t sumx3=0;
	Float_t sumx4=0;
	Float_t sumy=0;
	Float_t sumxy=0;
	Float_t sumx2y=0;
	Float_t sumz=0;
	Float_t sumxz=0;
	Float_t sumx2z=0;
	Float_t sumq=0;
	for (Int_t index=0;index<currentIndex;index++){
	  Float_t x,x2,x3,x4;
	  x=x2=x3=x4=xxx(index*4);
	  x2*=x;
	  x3*=x2;
	  x4*=x3;
	  sumx+=x;
	  sumx2+=x2;
	  sumx3+=x3;
	  sumx4+=x4;
	  sumy+=xxx(index*4+1);
	  sumxy+=xxx(index*4+1)*x;
	  sumx2y+=xxx(index*4+1)*x2;
	  sumz+=xxx(index*4+2);
	  sumxz+=xxx(index*4+2)*x;
	  sumx2z+=xxx(index*4+2)*x2;	 
	  sumq+=xxx(index*4+3);
	}
	Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
	Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
	  sumx2*(sumx*sumx3-sumx2*sumx2);
	
	Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
	  sumx2*(sumxy*sumx3-sumx2y*sumx2);
	Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
	  sumx2*(sumxz*sumx3-sumx2z*sumx2);
	
	Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
	  sumx2*(sumx*sumx2y-sumx2*sumxy);
	Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
	  sumx2*(sumx*sumx2z-sumx2*sumxz);
	
	Float_t y=detay/det+centralPad;
	Float_t z=detaz/det;	
	Float_t by=detby/det; //y angle
	Float_t bz=detbz/det; //z angle
	sumy/=Float_t(currentIndex);
	sumz/=Float_t(currentIndex);
	
	AliComplexCluster cl;
	cl.fX=z;
	cl.fY=y;
	cl.fQ=sumq;
	cl.fSigmaX2=bz;
	cl.fSigmaY2=by;
	cl.fTracks[0]=ipart;
	
	AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
	if (row!=0) row->InsertCluster(&cl);
	currentIndex=0;
	lastrow=currentrow;
      } //end of calculating cluster for given row
        	
  } // end of loop over hits
  xxxx->Delete();

}
Commit	Line	Data
6d75e4b6	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$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// fast TPC cluster simulation //
	21	// //
	22	///////////////////////////////////////////////////////////////////////////////
	23
	24	#include <TParticle.h>
e8d02863	25	#include <TVector.h>
6d75e4b6	26
	27	#include "AliRunLoader.h"
	28	#include "AliRun.h"
	29	#include "AliMC.h"
	30	#include "AliTPC.h"
	31	#include "AliTPCParam.h"
	32	#include "AliTPCClustersArray.h"
	33	#include "AliTPCClustersRow.h"
	34	#include "AliTPCcluster.h"
	35	#include "AliComplexCluster.h"
	36	#include "AliTPCFast.h"
8c2b3fd7	37	#include "AliLog.h"
6d75e4b6	38
	39	ClassImp(AliTPCFast)
	40
	41
	42	//_____________________________________________________________________________
	43	void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
	44	{
	45	//--------------------------------------------------------
	46	// TPC simple cluster generator from hits
	47	// obtained from the TPC Fast Simulator
	48	// The point errors are taken from the parametrization
	49	//--------------------------------------------------------
	50
	51	//-----------------------------------------------------------------
	52	// Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
	53	//-----------------------------------------------------------------
	54	// Adopted to Marian's cluster data structure by I.Belikov, CERN,
	55	// Jouri.Belikov@cern.ch
	56	//----------------------------------------------------------------
	57
	58	/////////////////////////////////////////////////////////////////////////////
	59	//
	60	//---------------------------------------------------------------------
	61	// ALICE TPC Cluster Parameters
	62	//--------------------------------------------------------------------
	63
	64
	65
	66	// Cluster width in rphi
	67	const Float_t kACrphi=0.18322;
	68	const Float_t kBCrphi=0.59551e-3;
	69	const Float_t kCCrphi=0.60952e-1;
	70	// Cluster width in z
	71	const Float_t kACz=0.19081;
	72	const Float_t kBCz=0.55938e-3;
	73	const Float_t kCCz=0.30428;
	74
	75
	76	AliLoader* loader = runLoader->GetLoader("TPCLoader");
	77	if (!loader) {
8c2b3fd7	78	AliError("No TPC loader found");
6d75e4b6	79	return;
	80	}
	81	if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
	82	AliRun* aliRun = runLoader->GetAliRun();
	83	if (!aliRun) {
8c2b3fd7	84	AliError("Couldn't get AliRun object");
6d75e4b6	85	return;
	86	}
	87	AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
	88	if (!tpc) {
8c2b3fd7	89	AliError("Couldn't get TPC detector");
6d75e4b6	90	return;
	91	}
	92	AliTPCParam* param = tpc->GetParam();
	93	if (!param) {
8c2b3fd7	94	AliError("No TPC parameters available");
6d75e4b6	95	return;
	96	}
	97
	98	//if(fDefaults == 0) SetDefaults();
	99
	100	Float_t sigmaRphi,sigmaZ,clRphi,clZ;
	101	//
	102	TParticle *particle; // pointer to a given particle
	103	AliTPChit *tpcHit; // pointer to a sigle TPC hit
	104	Int_t sector;
	105	Int_t ipart;
	106	Float_t xyz[5];
	107	Float_t pl,pt,tanth,rpad,ratio;
	108	Float_t cph,sph;
	109
	110	//---------------------------------------------------------------
	111	// Get the access to the tracks
	112	//---------------------------------------------------------------
	113
	114	TTree *tH = loader->TreeH();
	115	if (tH == 0x0)
8c2b3fd7	116	AliFatal("Can not find TreeH in folder");
8c2b3fd7	117
6d75e4b6	118	tpc->SetTreeAddress();
	119
	120	Stat_t ntracks = tH->GetEntries();
	121
	122	//Switch to the output file
	123
	124	if (loader->TreeR() == 0x0) loader->MakeTree("R");
	125
8c2b3fd7	126	AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
6d75e4b6	127
	128	runLoader->CdGAFile();
	129	//param->Write(param->GetTitle());
	130
	131	AliTPCClustersArray carray;
	132	carray.Setup(param);
	133	carray.SetClusterType("AliTPCcluster");
	134	carray.MakeTree(loader->TreeR());
	135
	136	Int_t nclusters=0; //cluster counter
	137
	138	//------------------------------------------------------------
	139	// Loop over all sectors (72 sectors for 20 deg
	140	// segmentation for both lower and upper sectors)
	141	// Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
	142	// Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
	143	//
	144	// First cluster for sector 0 starts at "0"
	145	//------------------------------------------------------------
	146
	147	for(Int_t isec=0;isec<param->GetNSector();isec++){
	148	//MI change
	149	param->AdjustCosSin(isec,cph,sph);
	150
	151	//------------------------------------------------------------
	152	// Loop over tracks
	153	//------------------------------------------------------------
	154
	155	for(Int_t track=0;track<ntracks;track++){
	156	tpc->ResetHits();
	157	tH->GetEvent(track);
	158	//
	159	// Get number of the TPC hits
	160	//
	161	tpcHit = (AliTPChit*)tpc->FirstHit(-1);
	162
	163	// Loop over hits
	164	//
	165	while(tpcHit){
	166
	167	if (tpcHit->fQ == 0.) {
	168	tpcHit = (AliTPChit*) tpc->NextHit();
	169	continue; //information about track (I.Belikov)
	170	}
	171	sector=tpcHit->fSector; // sector number
	172
	173	if(sector != isec){
	174	tpcHit = (AliTPChit*) tpc->NextHit();
	175	continue;
	176	}
	177	ipart=tpcHit->Track();
	178	particle=aliRun->GetMCApp()->Particle(ipart);
	179	pl=particle->Pz();
	180	pt=particle->Pt();
	181	if(pt < 1.e-9) pt=1.e-9;
	182	tanth=pl/pt;
	183	tanth = TMath::Abs(tanth);
	184	rpad=TMath::Sqrt(tpcHit->X()tpcHit->X() + tpcHit->Y()tpcHit->Y());
	185	ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
	186
	187	// space-point resolutions
	188
	189	sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
	190	sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
191
192	// cluster widths
193
194	clRphi=kACrphi-kBCrphirpadtanth+kCCrphiratioratio;
195	clZ=kACz-kBCzrpadtanth+kCCztanthtanth;
196
197	// temporary protection
198
199	if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
200	if(sigmaZ < 0.) sigmaZ=0.4e-3;
201	if(clRphi < 0.) clRphi=2.5e-3;
202	if(clZ < 0.) clZ=2.5e-5;
203
204	//
205
206	//
207	// smearing --> rotate to the 1 (13) or to the 25 (49) sector,
208	// then the inaccuracy in a X-Y plane is only along Y (pad row)!
209	//
210	Float_t xprim= tpcHit->X()cph + tpcHit->Y()sph;
211	Float_t yprim=-tpcHit->X()sph + tpcHit->Y()cph;
212	xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
213	Float_t alpha=(isec < param->GetNInnerSector()) ?
214	param->GetInnerAngle() : param->GetOuterAngle();
215	Float_t ymax=xprimTMath::Tan(0.5alpha);
216	if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
217	xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
218	if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
219	xyz[2]=sigmaRphi; // fSigmaY2
220	xyz[3]=sigmaZ; // fSigmaZ2
221	xyz[4]=tpcHit->fQ; // q
222
223	AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
224	if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
225
226	Int_t tracks[3]={tpcHit->Track(), -1, -1};
227	AliTPCcluster cluster(tracks,xyz);
228
229	clrow->InsertCluster(&cluster); nclusters++;
230
231	tpcHit = (AliTPChit*)tpc->NextHit();
232
233
234	} // end of loop over hits
235
236	} // end of loop over tracks
237
238	Int_t nrows=param->GetNRow(isec);
239	for (Int_t irow=0; irow<nrows; irow++) {
240	AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
241	if (!clrow) continue;
242	carray.StoreRow(isec,irow);
243	carray.ClearRow(isec,irow);
244	}
245
246	} // end of loop over sectors
247
248	// cerr<<"Number of made clusters : "<<nclusters<<" \n";
249	loader->WriteRecPoints("OVERWRITE");
250
251
252	} // end of function
253
254	//_________________________________________________________________
255	void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
256	AliTPCClustersArray* clustersArray,
257	Int_t isec) const
258	{
259	//--------------------------------------------------------
260	//calculate exact cross point of track and given pad row
261	//resulting values are expressed in "digit" coordinata
262	//--------------------------------------------------------
263
264	//-----------------------------------------------------------------
265	// Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
266	//-----------------------------------------------------------------
267	//
268	if (clustersArray==0){
269	return;
270	}
271	AliLoader* loader = runLoader->GetLoader("TPCLoader");
272	if (!loader) {
8c2b3fd7	273	AliError("No TPC loader found");
6d75e4b6	274	return;
	275	}
	276	if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
	277	AliRun* aliRun = runLoader->GetAliRun();
	278	if (!aliRun) {
8c2b3fd7	279	AliError("Couldn't get AliRun object");
6d75e4b6	280	return;
	281	}
	282	AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
	283	if (!tpc) {
8c2b3fd7	284	AliError("Couldn't get TPC detector");
6d75e4b6	285	return;
	286	}
	287	AliTPCParam* param = tpc->GetParam();
	288	if (!param) {
8c2b3fd7	289	AliError("No TPC parameters available");
6d75e4b6	290	return;
	291	}
	292	//
	293	TParticle *particle; // pointer to a given particle
	294	AliTPChit *tpcHit; // pointer to a sigle TPC hit
	295	// Int_t sector,nhits;
	296	Int_t ipart;
	297	const Int_t kcmaxhits=30000;
8c2b3fd7	298	TVector * xxxx = new TVector(kcmaxhits*4);
8c2b3fd7	299	TVector & xxx = *xxxx;
6d75e4b6	300	Int_t maxhits = kcmaxhits;
	301	//construct array for each padrow
	302	for (Int_t i=0; i<param->GetNRow(isec);i++)
	303	clustersArray->CreateRow(isec,i);
	304
	305	//---------------------------------------------------------------
	306	// Get the access to the tracks
	307	//---------------------------------------------------------------
	308
	309	TTree *tH = loader->TreeH();
	310	if (tH == 0x0)
8c2b3fd7	311	AliFatal("Can not find TreeH in folder");
8c2b3fd7	312
6d75e4b6	313	tpc->SetTreeAddress();
	314
	315	Stat_t ntracks = tH->GetEntries();
	316	Int_t npart = aliRun->GetMCApp()->GetNtrack();
	317	//MI change
	318	TBranch * branch=0;
	319	if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
	320	else branch = tH->GetBranch("TPC");
	321
	322	//------------------------------------------------------------
	323	// Loop over tracks
	324	//------------------------------------------------------------
	325
	326	for(Int_t track=0;track<ntracks;track++){
	327	Bool_t isInSector=kTRUE;
	328	tpc->ResetHits();
	329	isInSector = tpc->TrackInVolume(isec,track);
	330	if (!isInSector) continue;
	331	//MI change
	332	branch->GetEntry(track); // get next track
	333	//
	334	// Get number of the TPC hits and a pointer
	335	// to the particles
	336	// Loop over hits
	337	//
	338	Int_t currentIndex=0;
	339	Int_t lastrow=-1; //last writen row
	340
	341	//M.I. changes
	342
	343	tpcHit = (AliTPChit*)tpc->FirstHit(-1);
	344	while(tpcHit){
	345
	346	Int_t sector=tpcHit->fSector; // sector number
	347	if(sector != isec){
	348	tpcHit = (AliTPChit*) tpc->NextHit();
	349	continue;
	350	}
	351
	352	ipart=tpcHit->Track();
	353	if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
	354
	355	//find row number
	356
	357	Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
	358	Int_t index[3]={1,isec,0};
	359	Int_t currentrow = param->GetPadRow(x,index) ;
	360	if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
	361	if (lastrow<0) lastrow=currentrow;
	362	if (currentrow==lastrow){
	363	if ( currentIndex>=maxhits){
	364	maxhits+=kcmaxhits;
	365	xxx.ResizeTo(4*maxhits);
	366	}
	367	xxx(currentIndex*4)=x[0];
	368	xxx(currentIndex*4+1)=x[1];
	369	xxx(currentIndex*4+2)=x[2];
	370	xxx(currentIndex*4+3)=tpcHit->fQ;
	371	currentIndex++;
	372	}
	373	else
	374	if (currentIndex>2){
	375	Float_t sumx=0;
	376	Float_t sumx2=0;
377	Float_t sumx3=0;
378	Float_t sumx4=0;
379	Float_t sumy=0;
380	Float_t sumxy=0;
381	Float_t sumx2y=0;
382	Float_t sumz=0;
383	Float_t sumxz=0;
384	Float_t sumx2z=0;
385	Float_t sumq=0;
386	for (Int_t index=0;index<currentIndex;index++){
387	Float_t x,x2,x3,x4;
388	x=x2=x3=x4=xxx(index*4);
389	x2*=x;
390	x3*=x2;
391	x4*=x3;
392	sumx+=x;
393	sumx2+=x2;
394	sumx3+=x3;
395	sumx4+=x4;
396	sumy+=xxx(index*4+1);
397	sumxy+=xxx(index4+1)x;
398	sumx2y+=xxx(index4+1)x2;
399	sumz+=xxx(index*4+2);
400	sumxz+=xxx(index4+2)x;
401	sumx2z+=xxx(index4+2)x2;
402	sumq+=xxx(index*4+3);
403	}
404	Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
405	Float_t det=currentIndex(sumx2sumx4-sumx3sumx3)-sumx(sumxsumx4-sumx2sumx3)+
406	sumx2(sumxsumx3-sumx2*sumx2);
407
408	Float_t detay=sumy(sumx2sumx4-sumx3sumx3)-sumx(sumxysumx4-sumx2ysumx3)+
409	sumx2(sumxysumx3-sumx2y*sumx2);
410	Float_t detaz=sumz(sumx2sumx4-sumx3sumx3)-sumx(sumxzsumx4-sumx2zsumx3)+
411	sumx2(sumxzsumx3-sumx2z*sumx2);
412
413	Float_t detby=currentIndex(sumxysumx4-sumx2ysumx3)-sumy(sumxsumx4-sumx2sumx3)+
414	sumx2(sumxsumx2y-sumx2*sumxy);
415	Float_t detbz=currentIndex(sumxzsumx4-sumx2zsumx3)-sumz(sumxsumx4-sumx2sumx3)+
416	sumx2(sumxsumx2z-sumx2*sumxz);
417
418	if (TMath::Abs(det)<0.00001){
419	tpcHit = (AliTPChit*)tpc->NextHit();
420	continue;
421	}
422
423	Float_t y=detay/det+centralPad;
424	Float_t z=detaz/det;
425	Float_t by=detby/det; //y angle
426	Float_t bz=detbz/det; //z angle
427	sumy/=Float_t(currentIndex);
428	sumz/=Float_t(currentIndex);
429
430	AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
431	if (row!=0) {
432	AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
433	// AliComplexCluster cl;
434	cl->fX=z;
435	cl->fY=y;
436	cl->fQ=sumq;
437	cl->fSigmaX2=bz;
438	cl->fSigmaY2=by;
439	cl->fTracks[0]=ipart;
440	}
441	currentIndex=0;
442	lastrow=currentrow;
443	} //end of calculating cluster for given row
444
445
446	tpcHit = (AliTPChit*)tpc->NextHit();
447	} // end of loop over hits
448	} // end of loop over tracks
449	//write padrows to tree
450	for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
451	clustersArray->StoreRow(isec,ii);
452	clustersArray->ClearRow(isec,ii);
453	}
454	xxxx->Delete();
455
456	}
457
458
459	void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
460	AliTPCClustersArray* clustersArray) const
461	{
462
463	//
464	//fill clones array with intersection of current point with the
465	//middle of the row
466	AliLoader* loader = runLoader->GetLoader("TPCLoader");
467	if (!loader) {
8c2b3fd7	468	AliError("No TPC loader found");
6d75e4b6	469	return;
	470	}
	471	if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
	472	AliRun* aliRun = runLoader->GetAliRun();
	473	if (!aliRun) {
8c2b3fd7	474	AliError("Couldn't get AliRun object");
6d75e4b6	475	return;
	476	}
	477	AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
	478	if (!tpc) {
8c2b3fd7	479	AliError("Couldn't get TPC detector");
6d75e4b6	480	return;
	481	}
	482	AliTPCParam* param = tpc->GetParam();
	483	if (!param) {
8c2b3fd7	484	AliError("No TPC parameters available");
6d75e4b6	485	return;
	486	}
	487
	488	Int_t sector;
	489	Int_t ipart;
	490
	491	const Int_t kcmaxhits=30000;
8c2b3fd7	492	TVector * xxxx = new TVector(kcmaxhits*4);
8c2b3fd7	493	TVector & xxx = *xxxx;
6d75e4b6	494	Int_t maxhits = kcmaxhits;
	495
	496	//
	497	AliTPChit * tpcHit=0;
	498	tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
	499	Int_t currentIndex=0;
	500	Int_t lastrow=-1; //last writen row
	501
	502	while (tpcHit){
	503	if (tpcHit==0) continue;
	504	sector=tpcHit->fSector; // sector number
	505	ipart=tpcHit->Track();
	506
	507	//find row number
	508
	509	Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
	510	Int_t index[3]={1,sector,0};
	511	Int_t currentrow = param->GetPadRow(x,index) ;
	512	if (currentrow<0) continue;
	513	if (lastrow<0) lastrow=currentrow;
	514	if (currentrow==lastrow){
	515	if ( currentIndex>=maxhits){
	516	maxhits+=kcmaxhits;
	517	xxx.ResizeTo(4*maxhits);
	518	}
	519	xxx(currentIndex*4)=x[0];
	520	xxx(currentIndex*4+1)=x[1];
	521	xxx(currentIndex*4+2)=x[2];
	522	xxx(currentIndex*4+3)=tpcHit->fQ;
	523	currentIndex++;
	524	}
	525	else
	526	if (currentIndex>2){
	527	Float_t sumx=0;
	528	Float_t sumx2=0;
	529	Float_t sumx3=0;
	530	Float_t sumx4=0;
	531	Float_t sumy=0;
	532	Float_t sumxy=0;
	533	Float_t sumx2y=0;
	534	Float_t sumz=0;
	535	Float_t sumxz=0;
	536	Float_t sumx2z=0;
	537	Float_t sumq=0;
	538	for (Int_t index=0;index<currentIndex;index++){
	539	Float_t x,x2,x3,x4;
	540	x=x2=x3=x4=xxx(index*4);
	541	x2*=x;
	542	x3*=x2;
	543	x4*=x3;
	544	sumx+=x;
	545	sumx2+=x2;
	546	sumx3+=x3;
	547	sumx4+=x4;
	548	sumy+=xxx(index*4+1);
	549	sumxy+=xxx(index4+1)x;
	550	sumx2y+=xxx(index4+1)x2;
	551	sumz+=xxx(index*4+2);
	552	sumxz+=xxx(index4+2)x;
	553	sumx2z+=xxx(index4+2)x2;
	554	sumq+=xxx(index*4+3);
	555	}
	556	Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
	557	Float_t det=currentIndex(sumx2sumx4-sumx3sumx3)-sumx(sumxsumx4-sumx2sumx3)+
558	sumx2(sumxsumx3-sumx2*sumx2);
559
560	Float_t detay=sumy(sumx2sumx4-sumx3sumx3)-sumx(sumxysumx4-sumx2ysumx3)+
561	sumx2(sumxysumx3-sumx2y*sumx2);
562	Float_t detaz=sumz(sumx2sumx4-sumx3sumx3)-sumx(sumxzsumx4-sumx2zsumx3)+
563	sumx2(sumxzsumx3-sumx2z*sumx2);
564
565	Float_t detby=currentIndex(sumxysumx4-sumx2ysumx3)-sumy(sumxsumx4-sumx2sumx3)+
566	sumx2(sumxsumx2y-sumx2*sumxy);
567	Float_t detbz=currentIndex(sumxzsumx4-sumx2zsumx3)-sumz(sumxsumx4-sumx2sumx3)+
568	sumx2(sumxsumx2z-sumx2*sumxz);
569
570	Float_t y=detay/det+centralPad;
571	Float_t z=detaz/det;
572	Float_t by=detby/det; //y angle
573	Float_t bz=detbz/det; //z angle
574	sumy/=Float_t(currentIndex);
575	sumz/=Float_t(currentIndex);
576
577	AliComplexCluster cl;
578	cl.fX=z;
579	cl.fY=y;
580	cl.fQ=sumq;
581	cl.fSigmaX2=bz;
582	cl.fSigmaY2=by;
583	cl.fTracks[0]=ipart;
584
585	AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
586	if (row!=0) row->InsertCluster(&cl);
587	currentIndex=0;
588	lastrow=currentrow;
589	} //end of calculating cluster for given row
590
591	} // end of loop over hits
592	xxxx->Delete();
593
594	}
595