[u/mrichter/AliRoot.git] / TRD / AliTRDAnalysisTaskTP.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: AliTRDAnalysisTaskTP.cxx 42548 2010-07-27 08:10:51Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Track point maker for the alignment of TRD                            //
//                                                                        //
//  Author:                                                               //
//     Sebastian Huber (S.Huber@gsi.de)                                   // 
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <iostream>

#include <TROOT.h>
#include "TH1D.h"
#include "TH2D.h"
#include "TChain.h"
#include <TLinearFitter.h>

#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliInputEventHandler.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"
#include "AliAlignObjParams.h"
#include "AliTrackPointArray.h"
#include "AliESDfriend.h"
#include "AliESDtrack.h"
#include "AliESDfriendTrack.h"
#include "TFile.h"
#include "TTree.h"
#include "TSystem.h"
#include "TTimeStamp.h"
#include "TCollection.h"
#include "AliLog.h"
#include "AliGeomManager.h"

#include "AliTRDAnalysisTaskTP.h"

ClassImp(AliTRDAnalysisTaskTP)

AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP()
  :AliAnalysisTaskSE(),
  fArrHists(0x0),
  fArrTTree(0x0),
  fTree(NULL),
  fESD(0x0),
  fModpop(0x0),
  fBug(0x0),
  fNevents(0x0),
  fNtracks(0x0),
  fNAcceptedTracks(0),
  fArray(0x0),
  fFile(0x0)
{
  // default constructor

}
//____________________________________________________________
AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
  AliAnalysisTaskSE(name),
  fArrHists(0x0),
  fArrTTree(0x0),
  fTree(0x0),
  fESD(0x0),
  fModpop(0x0),
  fBug(0x0),
  fNevents(0),
  fNtracks(0),
  fNAcceptedTracks(0),
  fArray(0x0),
  fFile(0x0)
{
  // Constructor
  DefineOutput(1, TTree::Class());
  DefineOutput(2, TObjArray::Class());

}
//____________________________________________________________
AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP() {
  // destructor

}
//____________________________________________________________
void AliTRDAnalysisTaskTP::UserCreateOutputObjects() {

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

if (!fArrHists) fArrHists=new TObjArray;

  fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
  fModpop->SetXTitle("module nr");
  fModpop->SetYTitle("layer nr");
  fArrHists->Add(fModpop);

  OpenFile(1);
  fTree = new TTree("spTree", "Tree with track space point arrays");
  fTree->Branch("SP","AliTrackPointArray", &fArray);

}
//____________________________________________________________
void AliTRDAnalysisTaskTP::UserExec(Option_t *) {

  //AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
  fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
  if(!fESD){
  //cout << "ERROR: fESDs not available " << endl;
  return;
  }

  fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
  if (!fESDfriend) {
  //cout << "ERROR: fESDfriends not available " << endl;
  return;
  }

  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;
  const Float_t kMinPt          = 0;
  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;

  for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop

    AliESDtrack * track = fESD->GetTrack(itrack);
    fNtracks++;
    if (!track) continue;

    if (track->GetP() < kMinPt) continue;
    if (track->GetKinkIndex(0)!=0) continue;
    if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
    if (ptu) if (track->GetP() < kMinPtLow) continue;
    if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;

    AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
    if (!friendtrack){ 
    continue;
    }

    array = friendtrack->GetTrackPointArray();
    if (!array) continue;
    Int_t npoints = array->GetNPoints();
    if (tmpArray) delete tmpArray;
    tmpArray = new AliTrackPointArray(npoints);
    Int_t current = 0;
    int ntpc = 0;
    int ntrd = 0;
    for (Int_t ipoint=0; ipoint<npoints; ipoint++){  

      AliTrackPoint p;
      array->GetPoint(p, ipoint);

      UShort_t volid = array->GetVolumeID()[ipoint];
      Int_t iModule;
      AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
      if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;

      if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;

      Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
      if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;


      if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
	if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;

	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;

	fitter.Eval();
	fitterz.Eval();
	Double_t fity =fitter.GetParameter(0); 
	Double_t fitz =fitterz.GetParameter(0); 
	if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
	if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
	if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
	ntpc++;
      }

      if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){

	trdarray->AddPoint(ntrdarray++,&p);
	fModpop->Fill(iModule,layer);
	ntrd++;
      }
      tmpArray->AddPoint(current,&p);
      current++;
    }
    if (ntpc < 100) continue;
    if (ntrd < 4) continue;
    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++;
  PostData(1,fTree);
  PostData(2,fArrHists);
}
//____________________________________________________________
void AliTRDAnalysisTaskTP::Terminate(Option_t */*option*/) {
  //
  // Terminate
  //
}
Commit	Line	Data
3937e89f	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: AliTRDAnalysisTaskTP.cxx 42548 2010-07-27 08:10:51Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Track point maker for the alignment of TRD //
	21	// //
	22	// Author: //
	23	// Sebastian Huber (S.Huber@gsi.de) //
	24	// //
	25	////////////////////////////////////////////////////////////////////////////
	26
	27	#include <iostream>
	28
	29	#include <TROOT.h>
	30	#include "TH1D.h"
	31	#include "TH2D.h"
	32	#include "TChain.h"
	33	#include <TLinearFitter.h>
	34
	35	#include "AliVEvent.h"
	36	#include "AliESDEvent.h"
	37	#include "AliInputEventHandler.h"
	38	#include "AliESDInputHandler.h"
	39	#include "AliAnalysisTaskSE.h"
	40	#include "AliAnalysisManager.h"
	41	#include "AliAlignObjParams.h"
	42	#include "AliTrackPointArray.h"
	43	#include "AliESDfriend.h"
	44	#include "AliESDtrack.h"
	45	#include "AliESDfriendTrack.h"
	46	#include "TFile.h"
	47	#include "TTree.h"
	48	#include "TSystem.h"
	49	#include "TTimeStamp.h"
	50	#include "TCollection.h"
	51	#include "AliLog.h"
	52	#include "AliGeomManager.h"
	53
	54	#include "AliTRDAnalysisTaskTP.h"
	55
	56	ClassImp(AliTRDAnalysisTaskTP)
	57
	58	AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP()
	59	:AliAnalysisTaskSE(),
	60	fArrHists(0x0),
	61	fArrTTree(0x0),
	62	fTree(NULL),
	63	fESD(0x0),
	64	fModpop(0x0),
65	fBug(0x0),
66	fNevents(0x0),
67	fNtracks(0x0),
68	fNAcceptedTracks(0),
69	fArray(0x0),
70	fFile(0x0)
71	{
72	// default constructor
73
74	}
75	//____________________________________________________________
76	AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
77	AliAnalysisTaskSE(name),
78	fArrHists(0x0),
79	fArrTTree(0x0),
80	fTree(0x0),
81	fESD(0x0),
82	fModpop(0x0),
83	fBug(0x0),
84	fNevents(0),
85	fNtracks(0),
86	fNAcceptedTracks(0),
87	fArray(0x0),
88	fFile(0x0)
89	{
90	// Constructor
91	DefineOutput(1, TTree::Class());
92	DefineOutput(2, TObjArray::Class());
93
94	}
95	//____________________________________________________________
96	AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP() {
97	// destructor
98
99	}
100	//____________________________________________________________
101	void AliTRDAnalysisTaskTP::UserCreateOutputObjects() {
102
103	AliAlignObjParams alobj; // initialize align obj.
104	TString option = GetOption();
105
106	if (!fArrHists) fArrHists=new TObjArray;
107
108	fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
109	fModpop->SetXTitle("module nr");
110	fModpop->SetYTitle("layer nr");
111	fArrHists->Add(fModpop);
112
113	OpenFile(1);
114	fTree = new TTree("spTree", "Tree with track space point arrays");
115	fTree->Branch("SP","AliTrackPointArray", &fArray);
116
117	}
118	//____________________________________________________________
119	void AliTRDAnalysisTaskTP::UserExec(Option_t *) {
120
121	//AliESDEvent fESD = dynamic_cast<AliESDEvent >(fInputEvent);
122	fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
123	if(!fESD){
124	//cout << "ERROR: fESDs not available " << endl;
125	return;
126	}
127
128	fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
129	if (!fESDfriend) {
130	//cout << "ERROR: fESDfriends not available " << endl;
131	return;
132	}
133
134	TLinearFitter fitter(2, "pol1");
135	TLinearFitter fitterz(2, "pol1");
136
137	// track cuts
138	int tpc = 0; // require tpc
139	int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
140
141	const Float_t kMaxDelta = 1;
142	const Float_t kMinNcl = 60;
143	const Float_t kMinPtLow = 0.2;
144	const Float_t kMinNclLow = 100;
145	const Float_t kMinPt0 = 2;
146	const Float_t kMinPt = 0;
147	UInt_t status = AliESDtrack::kTRDrefit;
148	if (tpc) status \|= AliESDtrack::kTPCrefit;
149
150	const Float_t kMinRadius2 = 2*2;
151	const Float_t kMaxRadius2 = 400*400;
152	const Float_t kDeadSpace = 4;
153	const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
154	Int_t ntracks = fESD->GetNumberOfTracks();
155	const AliTrackPointArray *array=0;
156	AliTrackPointArray *tmpArray = 0;
157	// trdarray contains all trd points in this event, for duplication detection
158	AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
159	int ntrdarray = 0;
160
161	for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop
162
163	AliESDtrack * track = fESD->GetTrack(itrack);
164	fNtracks++;
165	if (!track) continue;
166
167	if (track->GetP() < kMinPt) continue;
168	if (track->GetKinkIndex(0)!=0) continue;
169	if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
170	if (ptu) if (track->GetP() < kMinPtLow) continue;
171	if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
172
173	AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
174	if (!friendtrack){
175	continue;
176	}
177
178	array = friendtrack->GetTrackPointArray();
179	if (!array) continue;
180	Int_t npoints = array->GetNPoints();
181	if (tmpArray) delete tmpArray;
182	tmpArray = new AliTrackPointArray(npoints);
183	Int_t current = 0;
184	int ntpc = 0;
185	int ntrd = 0;
186	for (Int_t ipoint=0; ipoint<npoints; ipoint++){
187
188	AliTrackPoint p;
189	array->GetPoint(p, ipoint);
190
191	UShort_t volid = array->GetVolumeID()[ipoint];
192	Int_t iModule;
193	AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
194	if ((layer < AliGeomManager::kFirstLayer) \|\| (layer >= AliGeomManager::kLastLayer)) continue;
195
196	if ((iModule >= AliGeomManager::LayerSize(layer)) \|\| (iModule < 0)) continue;
197
198	Float_t r2 = p.GetX()p.GetX()+p.GetY()p.GetY();
199	if ( r2<kMinRadius2 \|\| r2 > kMaxRadius2 ) continue;
200
201
202
203	if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
204	if (p.GetCov()[0]<0 \|\| p.GetCov()[3]<0 \|\| p.GetCov()[5]<0) continue;
205
206	AliTrackPoint& plocal = p.MasterToLocal();
207	Double_t ylocal = plocal.GetY();
208	Double_t zlocal = plocal.GetZ();
209	Double_t xlocal = plocal.GetX();
210	Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
211	Int_t nclose=0;
212	fitter.ClearPoints();
213	fitterz.ClearPoints();
214	for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
215	if (jpoint<0 \|\| jpoint>=npoints) continue;
216	if (ipoint==jpoint) continue;
217	UShort_t volidL = array->GetVolumeID()[jpoint];
218	if (volidL!=volid) continue;
219	AliTrackPoint pc;
220	array->GetPoint(pc, jpoint);
221	AliTrackPoint &pcl= pc.MasterToLocal();
222	Double_t dx = pcl.GetX()-xlocal;
223	fitter.AddPoint(&dx,pcl.GetY(),1);
224	fitterz.AddPoint(&dx,pcl.GetZ(),1);
225	nclose++;
226	}
227	if (nclose<6) continue;
228
229	fitter.Eval();
230	fitterz.Eval();
231	Double_t fity =fitter.GetParameter(0);
232	Double_t fitz =fitterz.GetParameter(0);
233	if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
234	if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
235	if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
236	ntpc++;
237	}
238
239	if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
240
241	trdarray->AddPoint(ntrdarray++,&p);
242	fModpop->Fill(iModule,layer);
243	ntrd++;
244	}
245	tmpArray->AddPoint(current,&p);
246	current++;
247	}
248	if (ntpc < 100) continue;
249	if (ntrd < 4) continue;
250	if (fArray) delete fArray;
251	fArray = new AliTrackPointArray(current);
252	for (Int_t ipoint=0; ipoint<current; ipoint++){
253	AliTrackPoint p;
254	tmpArray->GetPoint(p, ipoint);
255	fArray->AddPoint(ipoint,&p);
256	}
257	fNAcceptedTracks++;
258	fTree->Fill();
259	}
260	delete trdarray;
261	fNevents++;
262	PostData(1,fTree);
263	PostData(2,fArrHists);
264	}
265	//____________________________________________________________
266	void AliTRDAnalysisTaskTP::Terminate(Option_t /option*/) {
267	//
268	// Terminate
269	//
270	}
271
272