[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() 
{
  //
  // Create the output objects
  //

  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 *) 
{
  //
  // Exec function
  //

  //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	{
04e58504	72	//
	73	// Default constructor
	74	//
3937e89f	75
3937e89f	76	}
04e58504	77
3937e89f	78	//____________________________________________________________
	79	AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
	80	AliAnalysisTaskSE(name),
	81	fArrHists(0x0),
	82	fArrTTree(0x0),
	83	fTree(0x0),
	84	fESD(0x0),
	85	fModpop(0x0),
	86	fBug(0x0),
	87	fNevents(0),
	88	fNtracks(0),
	89	fNAcceptedTracks(0),
	90	fArray(0x0),
	91	fFile(0x0)
	92	{
04e58504	93	//
3937e89f	94	// Constructor
04e58504	95	//
04e58504	96
3937e89f	97	DefineOutput(1, TTree::Class());
	98	DefineOutput(2, TObjArray::Class());
	99
	100	}
04e58504	101
3937e89f	102	//____________________________________________________________
04e58504	103	AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP()
	104	{
	105	//
3937e89f	106	// destructor
04e58504	107	//
3937e89f	108
3937e89f	109	}
04e58504	110
3937e89f	111	//____________________________________________________________
04e58504	112	void AliTRDAnalysisTaskTP::UserCreateOutputObjects()
	113	{
	114	//
	115	// Create the output objects
	116	//
3937e89f	117
	118	AliAlignObjParams alobj; // initialize align obj.
	119	TString option = GetOption();
	120
04e58504	121	if (!fArrHists) fArrHists=new TObjArray;
3937e89f	122
	123	fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
	124	fModpop->SetXTitle("module nr");
	125	fModpop->SetYTitle("layer nr");
	126	fArrHists->Add(fModpop);
	127
	128	OpenFile(1);
	129	fTree = new TTree("spTree", "Tree with track space point arrays");
	130	fTree->Branch("SP","AliTrackPointArray", &fArray);
	131
	132	}
04e58504	133
3937e89f	134	//____________________________________________________________
04e58504	135	void AliTRDAnalysisTaskTP::UserExec(Option_t *)
	136	{
	137	//
	138	// Exec function
	139	//
3937e89f	140
	141	//AliESDEvent fESD = dynamic_cast<AliESDEvent >(fInputEvent);
	142	fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
	143	if(!fESD){
	144	//cout << "ERROR: fESDs not available " << endl;
	145	return;
	146	}
	147
	148	fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
	149	if (!fESDfriend) {
	150	//cout << "ERROR: fESDfriends not available " << endl;
	151	return;
	152	}
	153
	154	TLinearFitter fitter(2, "pol1");
	155	TLinearFitter fitterz(2, "pol1");
	156
	157	// track cuts
	158	int tpc = 0; // require tpc
	159	int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
	160
	161	const Float_t kMaxDelta = 1;
	162	const Float_t kMinNcl = 60;
	163	const Float_t kMinPtLow = 0.2;
	164	const Float_t kMinNclLow = 100;
	165	const Float_t kMinPt0 = 2;
	166	const Float_t kMinPt = 0;
	167	UInt_t status = AliESDtrack::kTRDrefit;
	168	if (tpc) status \|= AliESDtrack::kTPCrefit;
	169
	170	const Float_t kMinRadius2 = 2*2;
	171	const Float_t kMaxRadius2 = 400*400;
	172	const Float_t kDeadSpace = 4;
	173	const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
	174	Int_t ntracks = fESD->GetNumberOfTracks();
	175	const AliTrackPointArray *array=0;
	176	AliTrackPointArray *tmpArray = 0;
	177	// trdarray contains all trd points in this event, for duplication detection
	178	AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
	179	int ntrdarray = 0;
	180
	181	for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop
	182
	183	AliESDtrack * track = fESD->GetTrack(itrack);
	184	fNtracks++;
	185	if (!track) continue;
	186
	187	if (track->GetP() < kMinPt) continue;
	188	if (track->GetKinkIndex(0)!=0) continue;
	189	if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
	190	if (ptu) if (track->GetP() < kMinPtLow) continue;
	191	if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
	192
	193	AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
	194	if (!friendtrack){
	195	continue;
	196	}
	197
	198	array = friendtrack->GetTrackPointArray();
	199	if (!array) continue;
	200	Int_t npoints = array->GetNPoints();
	201	if (tmpArray) delete tmpArray;
	202	tmpArray = new AliTrackPointArray(npoints);
	203	Int_t current = 0;
204	int ntpc = 0;
205	int ntrd = 0;
206	for (Int_t ipoint=0; ipoint<npoints; ipoint++){
207
208	AliTrackPoint p;
209	array->GetPoint(p, ipoint);
210
211	UShort_t volid = array->GetVolumeID()[ipoint];
212	Int_t iModule;
213	AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
214	if ((layer < AliGeomManager::kFirstLayer) \|\| (layer >= AliGeomManager::kLastLayer)) continue;
215
216	if ((iModule >= AliGeomManager::LayerSize(layer)) \|\| (iModule < 0)) continue;
217
218	Float_t r2 = p.GetX()p.GetX()+p.GetY()p.GetY();
219	if ( r2<kMinRadius2 \|\| r2 > kMaxRadius2 ) continue;
220
221
222
223	if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
224	if (p.GetCov()[0]<0 \|\| p.GetCov()[3]<0 \|\| p.GetCov()[5]<0) continue;
225
226	AliTrackPoint& plocal = p.MasterToLocal();
227	Double_t ylocal = plocal.GetY();
228	Double_t zlocal = plocal.GetZ();
229	Double_t xlocal = plocal.GetX();
230	Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
231	Int_t nclose=0;
232	fitter.ClearPoints();
233	fitterz.ClearPoints();
234	for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
235	if (jpoint<0 \|\| jpoint>=npoints) continue;
236	if (ipoint==jpoint) continue;
237	UShort_t volidL = array->GetVolumeID()[jpoint];
238	if (volidL!=volid) continue;
239	AliTrackPoint pc;
240	array->GetPoint(pc, jpoint);
241	AliTrackPoint &pcl= pc.MasterToLocal();
242	Double_t dx = pcl.GetX()-xlocal;
243	fitter.AddPoint(&dx,pcl.GetY(),1);
244	fitterz.AddPoint(&dx,pcl.GetZ(),1);
245	nclose++;
246	}
247	if (nclose<6) continue;
248
249	fitter.Eval();
250	fitterz.Eval();
251	Double_t fity =fitter.GetParameter(0);
252	Double_t fitz =fitterz.GetParameter(0);
253	if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
254	if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
255	if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
256	ntpc++;
257	}
258
259	if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
260
261	trdarray->AddPoint(ntrdarray++,&p);
262	fModpop->Fill(iModule,layer);
263	ntrd++;
264	}
265	tmpArray->AddPoint(current,&p);
266	current++;
267	}
268	if (ntpc < 100) continue;
269	if (ntrd < 4) continue;
270	if (fArray) delete fArray;
271	fArray = new AliTrackPointArray(current);
272	for (Int_t ipoint=0; ipoint<current; ipoint++){
273	AliTrackPoint p;
274	tmpArray->GetPoint(p, ipoint);
275	fArray->AddPoint(ipoint,&p);
276	}
277	fNAcceptedTracks++;
278	fTree->Fill();
279	}
280	delete trdarray;
281	fNevents++;
282	PostData(1,fTree);
283	PostData(2,fArrHists);
284	}
04e58504	285
3937e89f	286	//____________________________________________________________
04e58504	287	void AliTRDAnalysisTaskTP::Terminate(Option_t /option*/)
04e58504	288	{
3937e89f	289	//
	290	// Terminate
	291	//
04e58504	292
3937e89f	293	}
	294
	295