[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDpointMaker.C

#define PointMaker_cxx
// The class definition in esdTree.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called everytime a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//

// To do: rewrite as AliAnalysisTaskSE.(like initiated by Alex at some point)

#include <TROOT.h>
#include <TStyle.h>
#include <TLinearFitter.h>
#include "AliAlignObjParams.h"
#include "AliTrackPointArray.h"
#include "AliLog.h"
#include "PointMaker.h"
//=============================================================================
PointMaker::PointMaker(char *outfil) :
  TSelector(),
  fChain(0),
  fESD(0),
  //  fESDfriend(0),
  fFile(0),
  fTree(0),
  fArray(0),
  fNevents(0),
  fNtracks(0),
  fNAcceptedTracks(0),
  fOutfil(outfil)
{
  // Constructor. Initialization of pointers
}
//=============================================================================
PointMaker::~PointMaker() {
  // Remove all pointers

  //  delete fESD;

  // histograms are in the output list and deleted when the output
  // list is deleted by the TSelector dtor
  delete fFile;
}
//=============================================================================
void PointMaker::Begin(TTree *)
{
  // The Begin() function is called at the start of the query.
  // When running with PROOF Begin() is only called on the client.
  // The tree argument is deprecated (on PROOF 0 is passed).

  TString option = GetOption();

  // create output file and tree  - TEMPORARY -until root trees are merged in the memory
  fFile = TFile::Open(fOutfil.Data(), "RECREATE",0);
  fTree = new TTree("spTree", "Tree with track space point arrays");
  fTree->Branch("SP","AliTrackPointArray", &fArray);
  printf("Begin called\n");
}
//=============================================================================
void PointMaker::SlaveBegin(TTree * tree)
{
  // The SlaveBegin() function is called after the Begin() function.
  // When running with PROOF SlaveBegin() is called on each slave server.
  // The tree argument is deprecated (on PROOF 0 is passed).

  Init(tree);

  AliAlignObjParams alobj;  // initialize align obj.  
  TString option = GetOption();

  // histograms to monitor cut efficiency and module population
  fCuttra = new TH1D("cuttra","cuttra",20,0.5,20.5); 
  fCutpoi = new TH1D("cutpoi","cutpoi",20,0.5,20.5);
  fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
  fModpop->SetXTitle("module nr");
  fModpop->SetYTitle("layer nr");
  printf("SlaveBegin called\n");
}
//=============================================================================
void PointMaker::Init(TTree *tree)
{
  // The Init() function is called when the selector needs to initialize
  // a new tree or chain. Typically here the branch addresses of the tree
  // will be set. It is normaly not necessary to make changes to the
  // generated code, but the routine can be extended by the user if needed.
  // Init() will be called many times when running with PROOF.

  // Set branch addresses
  if (tree == 0) return;
  fChain = tree;
  
  /*
  fChain->SetBranchAddress("ESD",&fESD);
  fChain->SetBranchStatus("ESDfriend*",1);
  fChain->SetBranchAddress("ESDfriend.",&fESDfriend);
  */
  
  fESD = new AliESDEvent();
  fChain->SetBranchStatus("ESDfriend*",1);
  fESD->ReadFromTree(fChain);
  fESDfriend = (AliESDfriend*)fESD->FindListObject("AliESDfriend");
  if(!fESDfriend) fChain->SetBranchAddress("ESDfriend.",&fESDfriend); 
  printf("Init called\n");
}
//=============================================================================
Bool_t PointMaker::Notify()
{
  // The Notify() function is called when a new file is opened. This
  // can be either for a new TTree in a TChain or when when a new TTree
  // is started when using PROOF. Typically here the branch pointers
  // will be retrieved. It is normaly not necessary to make changes
  // to the generated code, but the routine can be extended by the
  // user if needed.

  printf("Notify called\n");
  return kTRUE;
}
//=============================================================================
Bool_t PointMaker::IsIdenticalWithOneOf(AliTrackPoint *p, AliTrackPointArray *parray, int nmax) {

  // Is the point p identical with one of the points on the list parray?
  // This is a fix for aliroot 4-16-Rev-01 (and before) writing some 
  // spurious uninitialized points. 

  for (int i=0; i<parray->GetNPoints() && i<nmax; i++) {
    AliTrackPoint pa;
    parray->GetPoint(pa,i);
    //printf("comparing %7.3f with %7.3f\n",p->GetY(),pa.GetY());
    if (p->GetResidual(pa,0)<1e-8) return kTRUE;
    //printf("different\n");
  }
  return kFALSE;
}
//=============================================================================
Bool_t PointMaker::Process(Long64_t entry)
{
  // The Process() function is called for each entry in the tree (or possibly
  // keyed object in the case of PROOF) to be processed. The entry argument
  // specifies which entry in the currently loaded tree is to be processed.
  // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
  // to read either all or the required parts of the data. When processing
  // keyed objects with PROOF, the object is already loaded and is available
  // via the fObject pointer.
  //
  // This function should contain the "body" of the analysis. It can contain
  // simple or elaborate selection criteria, run algorithms on the data
  // of the event and typically fill histograms.

  // WARNING when a selector is used with a TChain, you must use
  //  the pointer to the current TTree to call GetEntry(entry).
  //  The entry is always the local entry number in the current tree.
  //  Assuming that fChain is the pointer to the TChain being processed,
  //  use fChain->GetTree()->GetEntry(entry).

  fChain->GetTree()->GetEntry(entry);
  fESD->SetESDfriend(fESDfriend);

  TLinearFitter fitter(2, "pol1");
  TLinearFitter fitterz(2, "pol1");

  // track cuts

  int tpc = 0; // require tpc
  int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)

  const Float_t kMaxDelta       = 1;
  const Float_t kMinNcl         = 60;
  const Float_t kMinPtLow       = 0.2;
  const Float_t kMinNclLow      = 100;
  const Float_t kMinPt0         = 2;
  UInt_t status = AliESDtrack::kTRDrefit; 
  if (tpc) status |= AliESDtrack::kTPCrefit; 

  const Float_t kMinRadius2  = 2*2;
  const Float_t kMaxRadius2  = 400*400;
  const Float_t kDeadSpace   = 4;
  const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
  Int_t ntracks = fESD->GetNumberOfTracks();
  const AliTrackPointArray *array = 0;
  AliTrackPointArray *tmpArray = 0;
  // trdarray contains all trd points in this event, for duplication detection
  AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
  int ntrdarray = 0;

  int pr = (!gROOT->IsBatch());
  for (Int_t itrack=0; itrack < ntracks; itrack++) {
    if (pr) printf("\revent %d   track %d",fNevents,itrack);
    fCuttra->Fill(1);
    AliESDtrack * track = fESD->GetTrack(itrack);
    fNtracks++;
    if (!track) continue;
    fCuttra->Fill(2);
    //if ((track->GetStatus() & status) == 0) continue;
    fCuttra->Fill(3);
    if (track->GetKinkIndex(0)!=0) continue;
    fCuttra->Fill(4);
    if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
    fCuttra->Fill(5);
    if (ptu) if (track->GetP() < kMinPtLow) continue;
    fCuttra->Fill(6);
    if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
    fCuttra->Fill(7);
    //
    // select points
    //
    array = track->GetTrackPointArray();    
    if (!array) continue;
    Int_t npoints = array->GetNPoints();
    if (tmpArray) delete tmpArray;
    tmpArray = new AliTrackPointArray(npoints);
    Int_t current = 0;
    int ntpc = 0; // number of good TPC points
    int ntrd = 0; // number of good TRD points
    for (Int_t ipoint=0; ipoint<npoints; ipoint++){
      fCutpoi->Fill(1);
      AliTrackPoint p;
      array->GetPoint(p, ipoint);

      // filter buggy points

      UShort_t volid = array->GetVolumeID()[ipoint];
      Int_t iModule;
      AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
      if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;
      fCutpoi->Fill(2);
      if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;
      fCutpoi->Fill(3);
      Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
      if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;
      fCutpoi->Fill(4);

      // TPC clusters handling

      if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
	// invalid covariance 
	if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;
	fCutpoi->Fill(7);

	// remove edge clusters

	AliTrackPoint& plocal = p.MasterToLocal();
	Double_t ylocal  = plocal.GetY();
	Double_t zlocal  = plocal.GetZ();
	Double_t xlocal  = plocal.GetX();
	Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
	Int_t nclose=0;
	fitter.ClearPoints();
	fitterz.ClearPoints();
	for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
	  if (jpoint<0 || jpoint>=npoints) continue;
	  if (ipoint==jpoint) continue;
	  UShort_t volidL = array->GetVolumeID()[jpoint];
	  if (volidL!=volid) continue;
	  AliTrackPoint pc;	
	  array->GetPoint(pc, jpoint);
	  AliTrackPoint &pcl=  pc.MasterToLocal();		  
	  Double_t dx = pcl.GetX()-xlocal;
	  fitter.AddPoint(&dx,pcl.GetY(),1);	  
	  fitterz.AddPoint(&dx,pcl.GetZ(),1);	  
	  nclose++;
	}
	if (nclose<6) continue;
	fCutpoi->Fill(8);
	fitter.Eval();
	fitterz.Eval();
	Double_t fity =fitter.GetParameter(0); 
	Double_t fitz =fitterz.GetParameter(0); 
	if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;  //outlier
	fCutpoi->Fill(9);
	if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;  //outlier
	fCutpoi->Fill(10);
	if (TMath::Abs(fity)>edgey-kDeadSpace) continue;   // remove edge clusters
	fCutpoi->Fill(11);
	ntpc++;
      }
      fCutpoi->Fill(12);

      // TRD track points handling

      if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
	fCutpoi->Fill(14);
	if (IsIdenticalWithOneOf(&p,trdarray,ntrdarray)) continue; // bug fix
	trdarray->AddPoint(ntrdarray++,&p);
	fCutpoi->Fill(15);
	fModpop->Fill(iModule,layer);
	ntrd++;
      }
      //printf("event %4d  track %4d  volid %4d  layer %4d  module %4d  %10.3f %10.3f %10.3f\n",
      //     (int)entry,itrack,volid,layer,iModule,p.GetX(),p.GetY(),atan2(p.GetY(),p.GetX())*180/3.1416);
      tmpArray->AddPoint(current,&p);
      current++;
      fCutpoi->Fill(16);
    }
    fCuttra->Fill(8);
    if (ntpc < 100) continue;
    fCuttra->Fill(9);
    if (ntrd < 4) continue;
    fCuttra->Fill(10);
    if (fArray) delete fArray;
    fArray = new AliTrackPointArray(current);
    for (Int_t ipoint=0; ipoint<current; ipoint++){
      AliTrackPoint p;
      tmpArray->GetPoint(p, ipoint);
      fArray->AddPoint(ipoint,&p);
    }   
    fNAcceptedTracks++;
    fTree->Fill();
  }
  delete trdarray;
  fNevents++;
  return kTRUE;
}
//=============================================================================
void PointMaker::SlaveTerminate()
{
  // The SlaveTerminate() function is called after all entries or objects
  // have been processed. When running with PROOF SlaveTerminate() is called
  // on each slave server.

  // Add the histograms to the output on each slave server
  AliInfo(Form("\nNumber of tracks:\tprocessed\t%d\taccepted\t%d", fNtracks, fNAcceptedTracks));
}
//=============================================================================
void PointMaker::Terminate()
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.

  //  fTree = dynamic_cast<TTree*>(fOutput->FindObject("spTree"));
  
  //  TFile* file = TFile::Open("AliTrackPoints.root", "RECREATE");
  fFile->cd();
  fTree->Write();
  fCuttra->Write();
  fCutpoi->Write();
  fModpop->Write();
  fModpop->Draw();
}
//=============================================================================
Commit	Line	Data
59fbd3a7	1	#define PointMaker_cxx
	2	// The class definition in esdTree.h has been generated automatically
	3	// by the ROOT utility TTree::MakeSelector(). This class is derived
	4	// from the ROOT class TSelector. For more information on the TSelector
	5	// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
	6
	7	// The following methods are defined in this file:
	8	// Begin(): called everytime a loop on the tree starts,
	9	// a convenient place to create your histograms.
	10	// SlaveBegin(): called after Begin(), when on PROOF called only on the
	11	// slave servers.
	12	// Process(): called for each event, in this function you decide what
	13	// to read and fill your histograms.
	14	// SlaveTerminate: called at the end of the loop on the tree, when on PROOF
	15	// called only on the slave servers.
	16	// Terminate(): called at the end of the loop on the tree,
	17	// a convenient place to draw/fit your histograms.
	18	//
	19
	20	// To do: rewrite as AliAnalysisTaskSE.(like initiated by Alex at some point)
	21
	22	#include <TROOT.h>
	23	#include <TStyle.h>
	24	#include <TLinearFitter.h>
	25	#include "AliAlignObjParams.h"
	26	#include "AliTrackPointArray.h"
	27	#include "AliLog.h"
	28	#include "PointMaker.h"
	29	//=============================================================================
	30	PointMaker::PointMaker(char *outfil) :
	31	TSelector(),
	32	fChain(0),
	33	fESD(0),
	34	// fESDfriend(0),
	35	fFile(0),
	36	fTree(0),
	37	fArray(0),
	38	fNevents(0),
	39	fNtracks(0),
	40	fNAcceptedTracks(0),
	41	fOutfil(outfil)
	42	{
	43	// Constructor. Initialization of pointers
	44	}
	45	//=============================================================================
	46	PointMaker::~PointMaker() {
	47	// Remove all pointers
	48
	49	// delete fESD;
	50
	51	// histograms are in the output list and deleted when the output
	52	// list is deleted by the TSelector dtor
	53	delete fFile;
	54	}
	55	//=============================================================================
	56	void PointMaker::Begin(TTree *)
	57	{
	58	// The Begin() function is called at the start of the query.
	59	// When running with PROOF Begin() is only called on the client.
	60	// The tree argument is deprecated (on PROOF 0 is passed).
	61
	62	TString option = GetOption();
	63
	64	// create output file and tree - TEMPORARY -until root trees are merged in the memory
65	fFile = TFile::Open(fOutfil.Data(), "RECREATE",0);
66	fTree = new TTree("spTree", "Tree with track space point arrays");
67	fTree->Branch("SP","AliTrackPointArray", &fArray);
68	printf("Begin called\n");
69	}
70	//=============================================================================
71	void PointMaker::SlaveBegin(TTree * tree)
72	{
73	// The SlaveBegin() function is called after the Begin() function.
74	// When running with PROOF SlaveBegin() is called on each slave server.
75	// The tree argument is deprecated (on PROOF 0 is passed).
76
77	Init(tree);
78
79	AliAlignObjParams alobj; // initialize align obj.
80	TString option = GetOption();
81
82	// histograms to monitor cut efficiency and module population
83	fCuttra = new TH1D("cuttra","cuttra",20,0.5,20.5);
84	fCutpoi = new TH1D("cutpoi","cutpoi",20,0.5,20.5);
85	fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
86	fModpop->SetXTitle("module nr");
87	fModpop->SetYTitle("layer nr");
88	printf("SlaveBegin called\n");
89	}
90	//=============================================================================
91	void PointMaker::Init(TTree *tree)
92	{
93	// The Init() function is called when the selector needs to initialize
94	// a new tree or chain. Typically here the branch addresses of the tree
95	// will be set. It is normaly not necessary to make changes to the
96	// generated code, but the routine can be extended by the user if needed.
97	// Init() will be called many times when running with PROOF.
98
99	// Set branch addresses
100	if (tree == 0) return;
101	fChain = tree;
102
103	/*
104	fChain->SetBranchAddress("ESD",&fESD);
105	fChain->SetBranchStatus("ESDfriend*",1);
106	fChain->SetBranchAddress("ESDfriend.",&fESDfriend);
107	*/
108
109	fESD = new AliESDEvent();
110	fChain->SetBranchStatus("ESDfriend*",1);
111	fESD->ReadFromTree(fChain);
112	fESDfriend = (AliESDfriend*)fESD->FindListObject("AliESDfriend");
113	if(!fESDfriend) fChain->SetBranchAddress("ESDfriend.",&fESDfriend);
114	printf("Init called\n");
115	}
116	//=============================================================================
117	Bool_t PointMaker::Notify()
118	{
119	// The Notify() function is called when a new file is opened. This
120	// can be either for a new TTree in a TChain or when when a new TTree
121	// is started when using PROOF. Typically here the branch pointers
122	// will be retrieved. It is normaly not necessary to make changes
123	// to the generated code, but the routine can be extended by the
124	// user if needed.
125
126	printf("Notify called\n");
127	return kTRUE;
128	}
129	//=============================================================================
130	Bool_t PointMaker::IsIdenticalWithOneOf(AliTrackPoint p, AliTrackPointArray parray, int nmax) {
131
132	// Is the point p identical with one of the points on the list parray?
133	// This is a fix for aliroot 4-16-Rev-01 (and before) writing some
134	// spurious uninitialized points.
135
136	for (int i=0; i<parray->GetNPoints() && i<nmax; i++) {
137	AliTrackPoint pa;
138	parray->GetPoint(pa,i);
139	//printf("comparing %7.3f with %7.3f\n",p->GetY(),pa.GetY());
140	if (p->GetResidual(pa,0)<1e-8) return kTRUE;
141	//printf("different\n");
142	}
143	return kFALSE;
144	}
145	//=============================================================================
146	Bool_t PointMaker::Process(Long64_t entry)
147	{
148	// The Process() function is called for each entry in the tree (or possibly
149	// keyed object in the case of PROOF) to be processed. The entry argument
150	// specifies which entry in the currently loaded tree is to be processed.
151	// It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
152	// to read either all or the required parts of the data. When processing
153	// keyed objects with PROOF, the object is already loaded and is available
154	// via the fObject pointer.
155	//
156	// This function should contain the "body" of the analysis. It can contain
157	// simple or elaborate selection criteria, run algorithms on the data
158	// of the event and typically fill histograms.
159
160	// WARNING when a selector is used with a TChain, you must use
161	// the pointer to the current TTree to call GetEntry(entry).
162	// The entry is always the local entry number in the current tree.
163	// Assuming that fChain is the pointer to the TChain being processed,
164	// use fChain->GetTree()->GetEntry(entry).
165
166	fChain->GetTree()->GetEntry(entry);
167	fESD->SetESDfriend(fESDfriend);
168
169	TLinearFitter fitter(2, "pol1");
170	TLinearFitter fitterz(2, "pol1");
171
172	// track cuts
173
174	int tpc = 0; // require tpc
175	int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
176
177	const Float_t kMaxDelta = 1;
178	const Float_t kMinNcl = 60;
179	const Float_t kMinPtLow = 0.2;
180	const Float_t kMinNclLow = 100;
181	const Float_t kMinPt0 = 2;
182	UInt_t status = AliESDtrack::kTRDrefit;
183	if (tpc) status \|= AliESDtrack::kTPCrefit;
184
185	const Float_t kMinRadius2 = 2*2;
186	const Float_t kMaxRadius2 = 400*400;
187	const Float_t kDeadSpace = 4;
188	const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
189	Int_t ntracks = fESD->GetNumberOfTracks();
190	const AliTrackPointArray *array = 0;
191	AliTrackPointArray *tmpArray = 0;
192	// trdarray contains all trd points in this event, for duplication detection
193	AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
194	int ntrdarray = 0;
195
196	int pr = (!gROOT->IsBatch());
197	for (Int_t itrack=0; itrack < ntracks; itrack++) {
198	if (pr) printf("\revent %d track %d",fNevents,itrack);
199	fCuttra->Fill(1);
200	AliESDtrack * track = fESD->GetTrack(itrack);
201	fNtracks++;
202	if (!track) continue;
203	fCuttra->Fill(2);
204	//if ((track->GetStatus() & status) == 0) continue;
205	fCuttra->Fill(3);
206	if (track->GetKinkIndex(0)!=0) continue;
207	fCuttra->Fill(4);
208	if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
209	fCuttra->Fill(5);
210	if (ptu) if (track->GetP() < kMinPtLow) continue;
211	fCuttra->Fill(6);
212	if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
213	fCuttra->Fill(7);
214	//
215	// select points
216	//
217	array = track->GetTrackPointArray();
218	if (!array) continue;
219	Int_t npoints = array->GetNPoints();
220	if (tmpArray) delete tmpArray;
221	tmpArray = new AliTrackPointArray(npoints);
222	Int_t current = 0;
223	int ntpc = 0; // number of good TPC points
224	int ntrd = 0; // number of good TRD points
225	for (Int_t ipoint=0; ipoint<npoints; ipoint++){
226	fCutpoi->Fill(1);
227	AliTrackPoint p;
228	array->GetPoint(p, ipoint);
229
230	// filter buggy points
231
232	UShort_t volid = array->GetVolumeID()[ipoint];
233	Int_t iModule;
234	AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
235	if ((layer < AliGeomManager::kFirstLayer) \|\| (layer >= AliGeomManager::kLastLayer)) continue;
236	fCutpoi->Fill(2);
237	if ((iModule >= AliGeomManager::LayerSize(layer)) \|\| (iModule < 0)) continue;
238	fCutpoi->Fill(3);
239	Float_t r2 = p.GetX()p.GetX()+p.GetY()p.GetY();
240	if ( r2<kMinRadius2 \|\| r2 > kMaxRadius2 ) continue;
241	fCutpoi->Fill(4);
242
243	// TPC clusters handling
244
245	if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
246	// invalid covariance
247	if (p.GetCov()[0]<0 \|\| p.GetCov()[3]<0 \|\| p.GetCov()[5]<0) continue;
248	fCutpoi->Fill(7);
249
250	// remove edge clusters
251
252	AliTrackPoint& plocal = p.MasterToLocal();
253	Double_t ylocal = plocal.GetY();
254	Double_t zlocal = plocal.GetZ();
255	Double_t xlocal = plocal.GetX();
256	Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
257	Int_t nclose=0;
258	fitter.ClearPoints();
259	fitterz.ClearPoints();
260	for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
261	if (jpoint<0 \|\| jpoint>=npoints) continue;
262	if (ipoint==jpoint) continue;
263	UShort_t volidL = array->GetVolumeID()[jpoint];
264	if (volidL!=volid) continue;
265	AliTrackPoint pc;
266	array->GetPoint(pc, jpoint);
267	AliTrackPoint &pcl= pc.MasterToLocal();
268	Double_t dx = pcl.GetX()-xlocal;
269	fitter.AddPoint(&dx,pcl.GetY(),1);
270	fitterz.AddPoint(&dx,pcl.GetZ(),1);
271	nclose++;
272	}
273	if (nclose<6) continue;
274	fCutpoi->Fill(8);
275	fitter.Eval();
276	fitterz.Eval();
277	Double_t fity =fitter.GetParameter(0);
278	Double_t fitz =fitterz.GetParameter(0);
279	if (TMath::Abs(ylocal-fity)>kMaxDelta) continue; //outlier
280	fCutpoi->Fill(9);
281	if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue; //outlier
282	fCutpoi->Fill(10);
283	if (TMath::Abs(fity)>edgey-kDeadSpace) continue; // remove edge clusters
284	fCutpoi->Fill(11);
285	ntpc++;
286	}
287	fCutpoi->Fill(12);
288
289	// TRD track points handling
290
291	if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
292	fCutpoi->Fill(14);
293	if (IsIdenticalWithOneOf(&p,trdarray,ntrdarray)) continue; // bug fix
294	trdarray->AddPoint(ntrdarray++,&p);
295	fCutpoi->Fill(15);
296	fModpop->Fill(iModule,layer);
297	ntrd++;
298	}
299	//printf("event %4d track %4d volid %4d layer %4d module %4d %10.3f %10.3f %10.3f\n",
300	// (int)entry,itrack,volid,layer,iModule,p.GetX(),p.GetY(),atan2(p.GetY(),p.GetX())*180/3.1416);
301	tmpArray->AddPoint(current,&p);
302	current++;
303	fCutpoi->Fill(16);
304	}
305	fCuttra->Fill(8);
306	if (ntpc < 100) continue;
307	fCuttra->Fill(9);
308	if (ntrd < 4) continue;
309	fCuttra->Fill(10);
310	if (fArray) delete fArray;
311	fArray = new AliTrackPointArray(current);
312	for (Int_t ipoint=0; ipoint<current; ipoint++){
313	AliTrackPoint p;
314	tmpArray->GetPoint(p, ipoint);
315	fArray->AddPoint(ipoint,&p);
316	}
317	fNAcceptedTracks++;
318	fTree->Fill();
319	}
320	delete trdarray;
321	fNevents++;
322	return kTRUE;
323	}
324	//=============================================================================
325	void PointMaker::SlaveTerminate()
326	{
327	// The SlaveTerminate() function is called after all entries or objects
328	// have been processed. When running with PROOF SlaveTerminate() is called
329	// on each slave server.
330
331	// Add the histograms to the output on each slave server
332	AliInfo(Form("\nNumber of tracks:\tprocessed\t%d\taccepted\t%d", fNtracks, fNAcceptedTracks));
333	}
334	//=============================================================================
335	void PointMaker::Terminate()
336	{
337	// The Terminate() function is the last function to be called during
338	// a query. It always runs on the client, it can be used to present
339	// the results graphically or save the results to file.
340
341	// fTree = dynamic_cast<TTree*>(fOutput->FindObject("spTree"));
342
343	// TFile* file = TFile::Open("AliTrackPoints.root", "RECREATE");
344	fFile->cd();
345	fTree->Write();
346	fCuttra->Write();
347	fCutpoi->Write();
348	fModpop->Write();
349	fModpop->Draw();
350	}
351	//=============================================================================