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

/// \class TestTPCTrackHits
/// Macro to compare TClonesArray hits with interpolated hits
/// ConvertHits1 read
///
/// \author MI

#include "alles.h"
#include "AliObjectArray.h"
#include "AliTPCTrackHits.h"
#include "AliArrayBranch.h"
#include "TestTPCTrackHits.h"
 
/// \cond CLASSIMP
ClassImp(AliTPChitD)
/// \endcond
 
void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits,  Bool_t debug, TClonesArray *arrd=0);
AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);

void ConvertHits(const char * benchmark, Bool_t debug)
{      

  /// convert - not parametrised hits stored in Clones array
  /// to

  TFile f("galice.root","update");  
  TClonesArray *arr = new TClonesArray("AliTPChit",100);
  TTree * treeCl =(TTree*)f.Get("TreeH0");
  TBranch *branch = treeCl->GetBranch("TPC"); 
  AliTPCTrackHits *  myhits = new AliTPCTrackHits;
  branch->SetAddress(&arr);
  //particle  
  TTree * treeP = (TTree*)f.Get("TreeK0");
  TBranch *branchP = treeP->GetBranch("Particles"); 
  TClonesArray *arrp = new TClonesArray("TParticle",100);
  branchP->SetAddress(&arrp);
  branchP->GetEvent(0);
  //
  TFile f2("treeh.root","recreate");
  f2.SetCompressionLevel(10);
  TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
  //treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
  treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits", 
	 &myhits,treeh2,4096,1));

  TFile f3("treehcl.root","recreate");
  f3.SetCompressionLevel(2);
  TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
  treeh3->Branch("TPC", &arr,64000,100);

  gBenchmark->Start(benchmark); 
  Int_t trsize = (Int_t)treeCl->GetEntries();
  for (Int_t i=0;i<300;i++){    
    arr->Clear(); 
    Int_t size = branch->GetEvent(i);     
    //if (size>0){
    if ((size>0)&& arr->GetEntriesFast()>0) {    
      f.cd();
      myhits->Clear();
      myhits =MakeTrack(arr,arrp,myhits);            
      CompareHits(arr,myhits,debug);
      f2.cd();
      treeh2->Fill();
      f3.cd();
      treeh3->Fill();      
      if ((i%10)==0) cout<<i<<"\n"; 
    }    
  }
  f3.cd();
  treeh3->Write();
  f2.cd();
  treeh2->Write();
  f2.Close();
  f.Close();
  delete myhits;
  gBenchmark->Stop(benchmark);
  gBenchmark->Show(benchmark); 
}  


void CompareHits(const char * benchmark, Bool_t debug)
{        
  /// compare persistent hits

  TFile f2("treeh.root");
  TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
  AliTPCTrackHits *  myhits = new AliTPCTrackHits ; 
  AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC"); 
  mybranch->SetAddress(&myhits);

    
  TClonesArray *arr = new TClonesArray("AliTPChit",100);
  TFile f3("treehcl.root");
  TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
  TBranch *branch = treeh3->GetBranch("TPC"); 
  branch->SetAddress(&arr);
   
  cout<<"Lets go!\n"; 
  gBenchmark->Start(benchmark);
  Int_t trsize = (Int_t)treeh3->GetEntries();
  for (Int_t i=0;i<300;i++){    
    Int_t size = branch->GetEvent(i);     
    mybranch->GetEvent(i);
    //if (arr){
    if ((arr->GetEntriesFast()>0) && size>0) {
      CompareHits(arr,myhits,debug);
      if ((i%10)==0) cout<<i<<"\n"; 
    }    
  }
  gBenchmark->Stop(benchmark);
  gBenchmark->Show(benchmark); 
}  


void CompareHitsG(const char * benchmark, Bool_t debug)
{        
  /// compare persistent hits

  TFile f2("galice.root");
  TTree * treeh2 = (TTree*)f2.Get("TreeH0");
  AliTPCTrackHits *  myhits = new AliTPCTrackHits ; 
  AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC2"); 
  mybranch->SetAddress(&myhits);

    
  TClonesArray *arr = new TClonesArray("AliTPChit",100);
  //TFile f3("treehcl.root");
  //TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
  TBranch *branch = treeh2->GetBranch("TPC"); 
  branch->SetAddress(&arr);

  TFile f3("treehdelta.root","recreate");
  f3.SetCompressionLevel(2);
  TTree * treeh3 = new TTree("DelataH","DeltaH");
  TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
  treeh3->Branch("TPC", &arrd,64000,100);
   
  cout<<"Lets go!\n"; 
  gBenchmark->Start(benchmark);
  Int_t trsize = treeh2->GetEntries();
  for (Int_t i=0;i<trsize;i++){    
    Int_t size = branch->GetEvent(i);     
    mybranch->GetEvent(i);
    //if (arr){
    if ((arr->GetEntriesFast()>0) && size>0) {
      CompareHits(arr,myhits,debug,arrd);      
    }  
    if ((i%10)==0) cout<<i<<"\n";
    treeh3->Fill();         
  }
  treeh3->Write();
  f3.Close();
  gBenchmark->Stop(benchmark);
  gBenchmark->Show(benchmark); 
}  


AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)  
{
  /// make track wit hits  according

  AliTPChit * hit;
  //  AliTPCTrackHits * myhits= new AliTPCTrackHits;
  myhits->SetHitPrecision(0.002);
  myhits->SetStepPrecision(0.003);  
  myhits->SetMaxDistance(100);  
  myhits->FlushHitStack(kTRUE);
  for (Int_t i=0;i<arr->GetEntriesFast();i++){
    hit = (AliTPChit*)arr->At(i);    
    if (hit){
      TParticle *p = (TParticle*)arrp->At(hit->GetTrack());
      Float_t momentum = TMath::Sqrt(p->Px()*p->Px()+p->Py()*p->Py());
      Float_t ran= 100.*gRandom->Rndm();
      if (ran<1.)  myhits->FlushHitStack(kTRUE);
      if (momentum<0.01) {//if small momentum - not precise
	myhits->SetHitPrecision(0.05);
	myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
			     hit->Z(), hit->fQ+1000);
      }
      else {
	myhits->SetHitPrecision(0.002);
	myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
			     hit->Z(), hit->fQ);
      }
    }
  }
  myhits->FlushHitStack();
  return myhits;  
}


void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
{     
  /// if debug option  kTRUE
  /// compare hits and write result  to the stdoutput

  AliTPChit * hit, *hit2;
  if (arrd) arrd->Clear();

  for (Int_t i=0;i<arr->GetEntriesFast();i++){
    hit = (AliTPChit*)arr->At(i);  
    if (hit){
      if (i==0) myhits->First();
      else myhits->Next();
      hit2 = myhits->GetHit();   
      if (!hit2) {
	hit2=0;
	if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
      }

      if (hit&&hit2){
	// if (hit){
	
	AliTrackHitsParam *param= myhits->GetParam();
	AliHitInfo * info = myhits->GetHitInfo();
	//	
	if (arrd) {
	  TClonesArray &larrd = *arrd;
	  AliTPChitD * h = new(larrd[i]) AliTPChitD;
	  h->SetX(hit->X());
	  h->SetY(hit->Y());
	  h->SetZ(hit->Z());
	  h->SetTrack(hit->GetTrack());
	  h->fQ = hit->fQ;
	  h->fSector = hit->fSector;
	  AliTPChit * hitd = h->GetDelta();
	  hitd->SetX(hit->X()-hit2->X());
	  hitd->SetY(hit->Y()-hit2->Y());
	  hitd->SetZ(hit->Z()-hit2->Z());
	  hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
	  hitd->fQ = hit->fQ-hit2->fQ;
	  hitd->fSector = hit->fSector-hit2->fSector;
	}

	if (debug){
	  Float_t dd = 
	    TMath::Sqrt(
			(hit->X()-hit2->X())*(hit->X()-hit2->X())+
			(hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
			(hit->Z()-hit2->Z())*(hit->Z()-hit2->Z())); 
	  printf("C1\t%d\t%d\t%d\t%d\t",
		 hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
	  printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
		 hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
	  printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
		 dd, param->fR,param->fZ,param->fFi,param->fAn,
		 param->fAd,param->fTheta);
	  printf("%d\t%d\t%d\n",
		 (Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
	}
      }        
    }
  }
} 


//Filter for paw visualisation of results
//
//sed filter  
 //sed '/*/d' dd.txt | sed '/C/d'|  cat -n > out.txt
 // sed -n '/C2/p' dd.txt | sed 's/C1//' | cat -n >out2.txt 
 // sed -n '/C3/p' dd.txt | sed 's/C2//' | cat -n >out3.txt 

//filter   
//myfilter C1 dd.txt >out1.txt
//myfilter C2 dd.txt >out2.txt
//myfilter C3 dd.txt >out3.txt
// sed -n 1,50000p
// sed -n 1,50000p dd.txt | myfilter C1  | sed -n 1,50000p >out1.txt

//myfilter C1 dd.txt | sed -n 1,50000p >out1.txt
//myfilter C2 dd.txt | sed -n 1,50000p >out2.txt
//myfilter C3 dd.txt | sed -n 1,50000p >out3.txt 
/*
paw visualisation
 Nt/Create 1 'Test of Match' 21 ! ! event v1 v2 t1 t2  x1 x2 y1  y2 z1 z2 dd r z fi an ad theta dist q1 q2
  Nt/Read 1 out1.txt 
 nt 
 Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r  z fi alfa theta dr1 ddr ddz ddrfi dd
  Nt/Read 2 out2.txt
 
 Nt/Create 3 'Test of Match2' 2 ! ! event stacksize 
  Nt/Read 3 out3.txt     
*/ 

void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
	    Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
	    Double_t fSumX4, Int_t n,
	    Double_t &a, Double_t &b, Double_t &c)
{
  /// recalc parameters not fixing origin point

  Double_t det = 
    n* (fSumX2*fSumX4-fSumX3*fSumX3) -
    fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
    fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
    
  if (TMath::Abs(det)>0){    
    a = 
      (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
       fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
       fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det; 
    b=
      (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
      fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
      fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
    c=
      (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
       fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
       fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;  
    cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
  }
}

void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
{
  Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
  fSumY = fSumYX = fSumYX2 = fSumX =  fSumX2 = fSumX3 =  fSumX4 =0;
  Double_t a2,b2,c2;
  for (Int_t i=0;i<n; i++){
    Float_t  x = gRandom->Rndm();
    Float_t  y = a+b*x+c*x*x;
    fSumY+=y;
    fSumYX+=y*x;
    fSumYX2+=y*x*x;
    fSumX += x;
    fSumX2 += x*x;
    fSumX3 += x*x*x;
    fSumX4 += x*x*x*x;
  }
  Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);  
}
Commit	Line	Data
b6895dd8	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
3a4edebe	16	/// \class TestTPCTrackHits
0fd1fc0b	17	/// Macro to compare TClonesArray hits with interpolated hits
3a4edebe	18	/// ConvertHits1 read
0fd1fc0b	19	///
0fd1fc0b	20	/// \author MI
cf98c13f	21
	22	#include "alles.h"
	23	#include "AliObjectArray.h"
	24	#include "AliTPCTrackHits.h"
	25	#include "AliArrayBranch.h"
	26	#include "TestTPCTrackHits.h"
	27
0fd1fc0b	28	/// \cond CLASSIMP
cf98c13f	29	ClassImp(AliTPChitD)
0fd1fc0b	30	/// \endcond
cf98c13f	31
	32	void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd=0);
	33	AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
	34
	35	void ConvertHits(const char * benchmark, Bool_t debug)
	36	{
cf98c13f	37
0fd1fc0b	38	/// convert - not parametrised hits stored in Clones array
0fd1fc0b	39	/// to
cf98c13f	40
	41	TFile f("galice.root","update");
	42	TClonesArray *arr = new TClonesArray("AliTPChit",100);
	43	TTree * treeCl =(TTree*)f.Get("TreeH0");
	44	TBranch *branch = treeCl->GetBranch("TPC");
	45	AliTPCTrackHits * myhits = new AliTPCTrackHits;
	46	branch->SetAddress(&arr);
	47	//particle
	48	TTree * treeP = (TTree*)f.Get("TreeK0");
	49	TBranch *branchP = treeP->GetBranch("Particles");
	50	TClonesArray *arrp = new TClonesArray("TParticle",100);
	51	branchP->SetAddress(&arrp);
	52	branchP->GetEvent(0);
	53	//
	54	TFile f2("treeh.root","recreate");
	55	f2.SetCompressionLevel(10);
	56	TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
	57	//treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
	58	treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits",
	59	&myhits,treeh2,4096,1));
	60
	61	TFile f3("treehcl.root","recreate");
	62	f3.SetCompressionLevel(2);
	63	TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
d0f40f23	64	treeh3->Branch("TPC", &arr,64000,100);
cf98c13f	65
	66	gBenchmark->Start(benchmark);
	67	Int_t trsize = (Int_t)treeCl->GetEntries();
	68	for (Int_t i=0;i<300;i++){
	69	arr->Clear();
	70	Int_t size = branch->GetEvent(i);
	71	//if (size>0){
	72	if ((size>0)&& arr->GetEntriesFast()>0) {
	73	f.cd();
	74	myhits->Clear();
	75	myhits =MakeTrack(arr,arrp,myhits);
	76	CompareHits(arr,myhits,debug);
	77	f2.cd();
	78	treeh2->Fill();
	79	f3.cd();
	80	treeh3->Fill();
	81	if ((i%10)==0) cout<<i<<"\n";
	82	}
	83	}
	84	f3.cd();
	85	treeh3->Write();
	86	f2.cd();
	87	treeh2->Write();
	88	f2.Close();
	89	f.Close();
	90	delete myhits;
	91	gBenchmark->Stop(benchmark);
	92	gBenchmark->Show(benchmark);
	93	}
	94
	95
	96
	97	void CompareHits(const char * benchmark, Bool_t debug)
	98	{
3a4edebe	99	/// compare persistent hits
3a4edebe	100
cf98c13f	101	TFile f2("treeh.root");
	102	TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
	103	AliTPCTrackHits * myhits = new AliTPCTrackHits ;
	104	AliObjectBranch mybranch = (AliObjectBranch)treeh2->GetBranch("TPC");
	105	mybranch->SetAddress(&myhits);
	106
	107
	108	TClonesArray *arr = new TClonesArray("AliTPChit",100);
	109	TFile f3("treehcl.root");
	110	TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
	111	TBranch *branch = treeh3->GetBranch("TPC");
	112	branch->SetAddress(&arr);
	113
	114	cout<<"Lets go!\n";
	115	gBenchmark->Start(benchmark);
	116	Int_t trsize = (Int_t)treeh3->GetEntries();
	117	for (Int_t i=0;i<300;i++){
	118	Int_t size = branch->GetEvent(i);
	119	mybranch->GetEvent(i);
	120	//if (arr){
	121	if ((arr->GetEntriesFast()>0) && size>0) {
	122	CompareHits(arr,myhits,debug);
	123	if ((i%10)==0) cout<<i<<"\n";
	124	}
	125	}
	126	gBenchmark->Stop(benchmark);
	127	gBenchmark->Show(benchmark);
	128	}
	129
	130
	131	void CompareHitsG(const char * benchmark, Bool_t debug)
	132	{
3a4edebe	133	/// compare persistent hits
3a4edebe	134
cf98c13f	135	TFile f2("galice.root");
	136	TTree * treeh2 = (TTree*)f2.Get("TreeH0");
	137	AliTPCTrackHits * myhits = new AliTPCTrackHits ;
	138	AliObjectBranch mybranch = (AliObjectBranch)treeh2->GetBranch("TPC2");
	139	mybranch->SetAddress(&myhits);
	140
	141
	142	TClonesArray *arr = new TClonesArray("AliTPChit",100);
	143	//TFile f3("treehcl.root");
	144	//TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
	145	TBranch *branch = treeh2->GetBranch("TPC");
	146	branch->SetAddress(&arr);
	147
	148	TFile f3("treehdelta.root","recreate");
	149	f3.SetCompressionLevel(2);
	150	TTree * treeh3 = new TTree("DelataH","DeltaH");
	151	TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
d0f40f23	152	treeh3->Branch("TPC", &arrd,64000,100);
cf98c13f	153
	154	cout<<"Lets go!\n";
	155	gBenchmark->Start(benchmark);
	156	Int_t trsize = treeh2->GetEntries();
	157	for (Int_t i=0;i<trsize;i++){
	158	Int_t size = branch->GetEvent(i);
	159	mybranch->GetEvent(i);
	160	//if (arr){
	161	if ((arr->GetEntriesFast()>0) && size>0) {
	162	CompareHits(arr,myhits,debug,arrd);
	163	}
	164	if ((i%10)==0) cout<<i<<"\n";
	165	treeh3->Fill();
	166	}
	167	treeh3->Write();
	168	f3.Close();
	169	gBenchmark->Stop(benchmark);
	170	gBenchmark->Show(benchmark);
	171	}
	172
	173
	174
	175
	176	AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)
	177	{
3a4edebe	178	/// make track wit hits according
3a4edebe	179
cf98c13f	180	AliTPChit * hit;
	181	// AliTPCTrackHits * myhits= new AliTPCTrackHits;
	182	myhits->SetHitPrecision(0.002);
	183	myhits->SetStepPrecision(0.003);
	184	myhits->SetMaxDistance(100);
	185	myhits->FlushHitStack(kTRUE);
	186	for (Int_t i=0;i<arr->GetEntriesFast();i++){
	187	hit = (AliTPChit*)arr->At(i);
	188	if (hit){
	189	TParticle p = (TParticle)arrp->At(hit->GetTrack());
	190	Float_t momentum = TMath::Sqrt(p->Px()p->Px()+p->Py()p->Py());
	191	Float_t ran= 100.*gRandom->Rndm();
	192	if (ran<1.) myhits->FlushHitStack(kTRUE);
	193	if (momentum<0.01) {//if small momentum - not precise
	194	myhits->SetHitPrecision(0.05);
	195	myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
	196	hit->Z(), hit->fQ+1000);
	197	}
	198	else {
	199	myhits->SetHitPrecision(0.002);
	200	myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
	201	hit->Z(), hit->fQ);
	202	}
	203	}
	204	}
	205	myhits->FlushHitStack();
	206	return myhits;
	207	}
	208
	209
	210
	211	void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
	212	{
3a4edebe	213	/// if debug option kTRUE
	214	/// compare hits and write result to the stdoutput
	215
cf98c13f	216	AliTPChit * hit, *hit2;
	217	if (arrd) arrd->Clear();
	218
	219	for (Int_t i=0;i<arr->GetEntriesFast();i++){
	220	hit = (AliTPChit*)arr->At(i);
	221	if (hit){
	222	if (i==0) myhits->First();
	223	else myhits->Next();
	224	hit2 = myhits->GetHit();
	225	if (!hit2) {
	226	hit2=0;
	227	if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
	228	}
	229
	230	if (hit&&hit2){
	231	// if (hit){
	232
	233	AliTrackHitsParam *param= myhits->GetParam();
	234	AliHitInfo * info = myhits->GetHitInfo();
	235	//
	236	if (arrd) {
	237	TClonesArray &larrd = *arrd;
	238	AliTPChitD * h = new(larrd[i]) AliTPChitD;
	239	h->SetX(hit->X());
	240	h->SetY(hit->Y());
	241	h->SetZ(hit->Z());
	242	h->SetTrack(hit->GetTrack());
	243	h->fQ = hit->fQ;
	244	h->fSector = hit->fSector;
	245	AliTPChit * hitd = h->GetDelta();
	246	hitd->SetX(hit->X()-hit2->X());
	247	hitd->SetY(hit->Y()-hit2->Y());
	248	hitd->SetZ(hit->Z()-hit2->Z());
	249	hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
	250	hitd->fQ = hit->fQ-hit2->fQ;
	251	hitd->fSector = hit->fSector-hit2->fSector;
	252	}
	253
	254	if (debug){
	255	Float_t dd =
	256	TMath::Sqrt(
	257	(hit->X()-hit2->X())*(hit->X()-hit2->X())+
	258	(hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
	259	(hit->Z()-hit2->Z())*(hit->Z()-hit2->Z()));
	260	printf("C1\t%d\t%d\t%d\t%d\t",
	261	hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
	262	printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
	263	hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
	264	printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
	265	dd, param->fR,param->fZ,param->fFi,param->fAn,
	266	param->fAd,param->fTheta);
	267	printf("%d\t%d\t%d\n",
	268	(Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
	269	}
	270	}
	271	}
	272	}
	273	}
	274
	275
	276
	277
	278	//Filter for paw visualisation of results
	279	//
280	//sed filter
281	//sed '/*/d' dd.txt \| sed '/C/d'\| cat -n > out.txt
282	// sed -n '/C2/p' dd.txt \| sed 's/C1//' \| cat -n >out2.txt
283	// sed -n '/C3/p' dd.txt \| sed 's/C2//' \| cat -n >out3.txt
284
285	//filter
286	//myfilter C1 dd.txt >out1.txt
287	//myfilter C2 dd.txt >out2.txt
288	//myfilter C3 dd.txt >out3.txt
289	// sed -n 1,50000p
290	// sed -n 1,50000p dd.txt \| myfilter C1 \| sed -n 1,50000p >out1.txt
291
292	//myfilter C1 dd.txt \| sed -n 1,50000p >out1.txt
293	//myfilter C2 dd.txt \| sed -n 1,50000p >out2.txt
294	//myfilter C3 dd.txt \| sed -n 1,50000p >out3.txt
295	/*
296	paw visualisation
297	Nt/Create 1 'Test of Match' 21 ! ! event v1 v2 t1 t2 x1 x2 y1 y2 z1 z2 dd r z fi an ad theta dist q1 q2
298	Nt/Read 1 out1.txt
299	nt
300	Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r z fi alfa theta dr1 ddr ddz ddrfi dd
301	Nt/Read 2 out2.txt
302
303	Nt/Create 3 'Test of Match2' 2 ! ! event stacksize
304	Nt/Read 3 out3.txt
305	*/
306
307	void Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
308	Double_t fSumX, Double_t fSumX2, Double_t fSumX3,
309	Double_t fSumX4, Int_t n,
310	Double_t &a, Double_t &b, Double_t &c)
311	{
3a4edebe	312	/// recalc parameters not fixing origin point
3a4edebe	313
cf98c13f	314	Double_t det =
	315	n* (fSumX2fSumX4-fSumX3fSumX3) -
	316	fSumX* (fSumXfSumX4-fSumX3fSumX2)+
	317	fSumX2* (fSumXfSumX3-fSumX2fSumX2);
	318
	319	if (TMath::Abs(det)>0){
	320	a =
	321	(fSumY * (fSumX2fSumX4-fSumX3fSumX3)-
	322	fSumX (fSumYXfSumX4-fSumYX2*fSumX3)+
	323	fSumX2(fSumYXfSumX3-fSumYX2*fSumX2))/det;
	324	b=
	325	(n(fSumYXfSumX4-fSumX3*fSumYX2)-
	326	fSumY(fSumXfSumX4-fSumX3*fSumX2)+
	327	fSumX2(fSumXfSumYX2-fSumYX*fSumX2))/det;
	328	c=
	329	(n(fSumX2fSumYX2-fSumYX*fSumX3)-
	330	fSumX(fSumXfSumYX2-fSumYX*fSumX2)+
	331	fSumY(fSumXfSumX3-fSumX2*fSumX2))/det;
	332	cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
	333	}
	334	}
	335
	336	void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
	337	{
	338	Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
	339	fSumY = fSumYX = fSumYX2 = fSumX = fSumX2 = fSumX3 = fSumX4 =0;
	340	Double_t a2,b2,c2;
	341	for (Int_t i=0;i<n; i++){
	342	Float_t x = gRandom->Rndm();
	343	Float_t y = a+bx+cx*x;
	344	fSumY+=y;
	345	fSumYX+=y*x;
	346	fSumYX2+=yxx;
	347	fSumX += x;
	348	fSumX2 += x*x;
	349	fSumX3 += xxx;
	350	fSumX4 += xxx*x;
	351	}
	352	Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);
	353	}