[u/mrichter/AliRoot.git] / HLT / MUON / macros / MakeTrackTable.C

/**************************************************************************
 * This file is property of and copyright by the ALICE HLT Project        * 
 * All rights reserved.                                                   *
 *                                                                        *
 * Primary Authors:                                                       *
 *   Artur Szostak <artursz@iafrica.com>                                  *
 *                                                                        *
 * 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.                  *
 **************************************************************************/

// Simple macro to generate N-tuple for performance analysis.
// This macro must be compiled and run like so from within the aliroot command
// prompt:
//   root [0] gSystem->Load("libAliHLTMUON");
//   root [0] .x MakeHitsTable.C+

#include "TVector3.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TParticle.h"
#include "TArrayF.h"
#include "TStopwatch.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TError.h"
#include "TMath.h"

#include "AliCDBManager.h"
#include "AliLoader.h"
#include "AliRunLoader.h"
#include "AliStack.h"
#include "AliRun.h"
#include "AliMUON.h"
#include "AliMUONMCDataInterface.h"
#include "AliMUONHit.h"
#include "AliHLTMUONEvent.h"
#include "AliHLTMUONRootifierComponent.h"
#include "AliHLTMUONMansoTrack.h"

#include <cstdlib>
#include <vector>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;


void MakeTrackTable(
		Int_t firstEvent = 0,
		Int_t lastEvent = -1,
		const char* dHLToutputfile = "output.root",
		Float_t maxSigma = 4., // 4 standard deviations
		Float_t sigmaX = 0.1,  // 1 mm resolution
		Float_t sigmaY = 0.01, // 100 micron resolution
		Float_t sigmaZ = 0.02,  // 200 microns resolution
		Float_t sigmaXtrg = 0.5,  // 5 mm resolution
		Float_t sigmaYtrg = 0.5,  // 5 mm resolution
		Float_t sigmaZtrg = 0.02  // 2 microns resolution
	)
{
	// Setup the CDB default storage and run number if nothing was set.
	AliCDBManager* cdbManager = AliCDBManager::Instance();
	if (cdbManager == NULL)
	{
		cerr << "ERROR: Global CDB manager object does not exist." << endl;
		return;
	}
	if (cdbManager->GetDefaultStorage() == NULL)
	{
		cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
	}
	if (cdbManager->GetRun() == -1)
	{
		cdbManager->SetRun(0);
	}
	
	gSystem->Load("libAliHLTMUON");
	// Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.

	TString fieldnames = "event:isprimary:hastrack:cantrigger:pdgcode:sign:px:py:pz";
	for (Int_t i = 0; i <= 14; i++)
	{
		fieldnames += ":x";
		fieldnames += i;
		fieldnames += ":y";
		fieldnames += i;
		fieldnames += ":z";
		fieldnames += i;
	}
	fieldnames += ":dHLTsign:dHLTpx:dHLTpy:dHLTpz:dHLTfit";
	for (Int_t i = 1; i <= 14; i++)
	{
		fieldnames += ":dHLTx";
		fieldnames += i;
		fieldnames += ":dHLTy";
		fieldnames += i;
		fieldnames += ":dHLTz";
		fieldnames += i;
	}
	TNtuple* table = new TNtuple("tracktable", "Table of tracks.", fieldnames);
	
	TFile dHLTfile(dHLToutputfile, "READ");

	AliMUONMCDataInterface data;
	AliStack* stack;
	
	Int_t nEvents = data.NumberOfEvents();
	if (lastEvent < 0) lastEvent = nEvents-1;
	for (Int_t event = firstEvent; event <= lastEvent; event++)
	{
		cout << "Processing event: " << event;
		
		char buf[1024];
		char* str = &buf[0];
		sprintf(str, "AliHLTMUONEvent;%d", event+1);
		AliHLTMUONEvent* dHLTevent = dynamic_cast<AliHLTMUONEvent*>( dHLTfile.Get(str) );
		if (dHLTevent == NULL)
		{
			cout << endl;
			cerr << "ERROR: Could not find " << str << " in " << dHLToutputfile << endl;
			continue;
		}
		
		data.GetEvent(event);
		
		stack = data.Stack(event);
		if (stack == NULL)
		{
			cout << endl;
			cerr << "ERROR: Stack is NULL" << endl;
			continue;
		}
		
		Int_t trackcount = data.NumberOfParticles();
		cout << "\ttrackcount = " << trackcount << endl;
		
		TArrayF isPrimary(trackcount);
		TArrayF hasTrack(trackcount);
		TArrayI hitsInTrigger(trackcount);
		TArrayF pdgcode(trackcount);
		TArrayF sign(trackcount);
		TArrayF px(trackcount);
		TArrayF py(trackcount);
		TArrayF pz(trackcount);
		TArrayF vx(trackcount);
		TArrayF vy(trackcount);
		TArrayF vz(trackcount);
		TArrayF dHLTsign(trackcount);
		TArrayF dHLTpx(trackcount);
		TArrayF dHLTpy(trackcount);
		TArrayF dHLTpz(trackcount);
		TArrayF dHLTfit(trackcount);
		
		//cout << "Fill particle data" << endl;
		for (Int_t i = 0; i < trackcount; i++)
		{
			TParticle* p = data.Particle(i);
			//isPrimary[i] = (i < stack->GetNprimary()) ? 1 : 0;
			isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
			hasTrack[i] = 0;
			hitsInTrigger[i] = 0;
			pdgcode[i] = p->GetPdgCode();
			if (p->GetPDG() != NULL)
			{
				Double_t charge = p->GetPDG()->Charge();
				if (charge > 0)
					sign[i] = 1;
				else if (charge < 0)
					sign[i] = -1;
				else
					sign[i] = 0;
			}
			else
			{
				sign[i] = 0;
				cerr << "ERROR: Could not get the PDG information." << endl;
			}
			px[i] = p->Px();
			py[i] = p->Py();
			pz[i] = p->Pz();
			vx[i] = p->Vx();
			vy[i] = p->Vy();
			vz[i] = p->Vz();
			dHLTsign[i] = 0;
			dHLTpx[i] = 0;
			dHLTpy[i] = 0;
			dHLTpz[i] = 0;
			dHLTfit[i] = -1;
		}
		
		TArrayF hitX[14];
		TArrayF hitY[14];
		TArrayF hitZ[14];
		TArrayF dHLThitX[14];
		TArrayF dHLThitY[14];
		TArrayF dHLThitZ[14];
		for (Int_t i = 0; i < 14; i++)
		{
			hitX[i].Set(trackcount);
			hitY[i].Set(trackcount);
			hitZ[i].Set(trackcount);
			dHLThitX[i].Set(trackcount);
			dHLThitY[i].Set(trackcount);
			dHLThitZ[i].Set(trackcount);
			for (Int_t j = 0; j < trackcount; j++)
			{
				hitX[i][j] = 0;
				hitY[i][j] = 0;
				hitZ[i][j] = 0;
				dHLThitX[i][j] = 0;
				dHLThitY[i][j] = 0;
				dHLThitZ[i][j] = 0;
			}
		}

		//cout << "Fill hits" << endl;
		for (Int_t i = 0; i < data.NumberOfTracks(); i++)
		for (Int_t j = 0; j < data.NumberOfHits(i); j++)
		{
			AliMUONHit* hit = data.Hit(i, j);
			Int_t chamber = hit->Chamber() - 1;
			if (chamber < 0 || chamber >= 14)
			{
				cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
				continue;
			}
			
			//cout << "hit->Track() = " << hit->Track() << endl;
			if (hit->Track() < 0 || hit->Track() >= trackcount)
			{
				cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
					<< trackcount << ")" << endl;
				continue;
			}
			Int_t particleindex = hit->Track();

			hitX[chamber][particleindex] = hit->X();
			hitY[chamber][particleindex] = hit->Y();
			hitZ[chamber][particleindex] = hit->Z();
					
			if (chamber == 0)
			{
				hasTrack[particleindex] = 1;
			}
			if (chamber >= 10)
			{
				hitsInTrigger[particleindex]++;
			}
		} // hits loop
		
		// Ok now go through the dHLT information for this event and correlate
		// the dHLT tracks to the kine data.
		// This is done as follows: take a dHLT track. Then for each kine track calculate
		// the distance between each hit of the kine track and corresponding dHLT track.
		// Then sum up the distances and divide by the number of hits matched.
		// This gives us a fit quality /Chi2 type parameter. The kine track that
		// has the smallest sum of differences is then chosen as the kine track
		// correlating to the current dHLT track being processed.
		for (Int_t i = 0; i < dHLTevent->Array().GetEntriesFast(); i++)
		{
			if (dHLTevent->Array().At(i)->IsA() != AliHLTMUONMansoTrack::Class())
				continue;
			AliHLTMUONMansoTrack* mtrack = static_cast<AliHLTMUONMansoTrack*>(dHLTevent->Array().At(i));
			
			// Find the best fitting kine track.
			Double_t bestFitQuality = 1e30;
			Int_t bestFitTrack = -1;
			for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++)
			{
				Double_t sumOfDiffs = 0;
				Double_t hitsMatched = 0;
				for (Int_t j = 7; j <= 10; j++)
				{
					const AliHLTMUONRecHit* hit = mtrack->Hit(j);
					if (hit == NULL) continue;
					if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
						continue;
					
					TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
					TVector3 diff = hV - hit->Coordinate();
					Double_t diffX = diff.X() / sigmaX;
					Double_t diffY = diff.Y() / sigmaY;
					Double_t diffZ = diff.Z() / sigmaZ;
					if (diffX > maxSigma) continue;
					if (diffY > maxSigma) continue;
					if (diffZ > maxSigma) continue;
					sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
					hitsMatched++;
				}
				if (mtrack->TriggerRecord() != NULL)
				{
					for (Int_t j = 11; j <= 14; j++)
					{
						const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
						if (hit == TVector3(0,0,0)) continue;
						if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
							continue;
						
						TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
						TVector3 diff = hV - hit;
						Double_t diffX = diff.X() / sigmaXtrg;
						Double_t diffY = diff.Y() / sigmaYtrg;
						Double_t diffZ = diff.Z() / sigmaZtrg;
						if (diffX > maxSigma) continue;
						if (diffY > maxSigma) continue;
						if (diffZ > maxSigma) continue;
						sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
						hitsMatched++;
					}
				}
				
				// Now check the fit quality.
				if (hitsMatched <= 0) continue;
				Double_t fitQuality = TMath::Sqrt(sumOfDiffs) / hitsMatched;
				if (fitQuality < bestFitQuality)
				{
					bestFitQuality = fitQuality;
					bestFitTrack = kinetrack;
				}
			}
			
			if (bestFitTrack < 0)
			{
				// No track was matched so we need to add a fake track record.
				Float_t args[12+14*3+5+14*3];
				args[0] = event;
				args[1] = -1;
				args[2] = 0;
				args[3] = 0;
				args[4] = 0;
				args[5] = 0;
				args[6] = 0;
				args[7] = 0;
				args[8] = 0;
				args[9] = 0;
				args[10] = 0;
				args[11] = 0;
				for (Int_t j = 0; j < 14; j++)
				{			
					args[12+j*3+0] = 0;
					args[12+j*3+1] = 0;
					args[12+j*3+2] = 0;
				}
				args[12+14*3+0] = mtrack->Sign();
				args[12+14*3+1] = mtrack->Px();
				args[12+14*3+2] = mtrack->Py();
				args[12+14*3+3] = mtrack->Pz();
				args[12+14*3+4] = -1;
				for (Int_t j = 0; j < 14; j++)
				{			
					args[12+14*3+5+j*3+0] = 0;
					args[12+14*3+5+j*3+1] = 0;
					args[12+14*3+5+j*3+2] = 0;
				}
				for (Int_t j = 6; j < 10; j++)
				{
					const AliHLTMUONRecHit* hit = mtrack->Hit(j+1);
					if (hit == NULL) continue;
					args[12+14*3+5+j*3+0] = hit->X();
					args[12+14*3+5+j*3+1] = hit->Y();
					args[12+14*3+5+j*3+2] = hit->Z();
				}
				if (mtrack->TriggerRecord() != NULL)
				{
					for (Int_t j = 10; j < 14; j++)
					{
						const TVector3& hit = mtrack->TriggerRecord()->Hit(j+1);
						args[12+14*3+5+j*3+0] = hit.X();
						args[12+14*3+5+j*3+1] = hit.Y();
						args[12+14*3+5+j*3+2] = hit.Z();
					}
				}
				table->Fill(args);
			}
			else
			{
				// Fill the details about the dHLT track to the best fitting track info.
				dHLTsign[bestFitTrack] = mtrack->Sign();
				dHLTpx[bestFitTrack] = mtrack->Px();
				dHLTpy[bestFitTrack] = mtrack->Py();
				dHLTpz[bestFitTrack] = mtrack->Pz();
				dHLTfit[bestFitTrack] = bestFitQuality;
				for (Int_t j = 7; j <= 10; j++)
				{
					const AliHLTMUONRecHit* hit = mtrack->Hit(j);
					if (hit == NULL) continue;
					dHLThitX[j-1][bestFitTrack] = hit->X();
					dHLThitY[j-1][bestFitTrack] = hit->Y();
					dHLThitZ[j-1][bestFitTrack] = hit->Z();
				}
				if (mtrack->TriggerRecord() != NULL)
				{
					for (Int_t j = 11; j <= 14; j++)
					{
						const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
						if (hit == TVector3(0,0,0)) continue;
						dHLThitX[j-1][bestFitTrack] = hit.X();
						dHLThitY[j-1][bestFitTrack] = hit.Y();
						dHLThitZ[j-1][bestFitTrack] = hit.Z();
					}
				}
			}
		}
		
		//cout << "Fill table" << endl;
		for (Int_t i = 0; i < trackcount; i++)
		{
			Float_t args[12+14*3+5+14*3];
			args[0] = event;
			args[1] = isPrimary[i];
			args[2] = hasTrack[i];
			args[3] = (hitsInTrigger[i] > 2) ? 1 : 0;
			args[4] = pdgcode[i];
			args[5] = sign[i];
			args[6] = px[i];
			args[7] = py[i];
			args[8] = pz[i];
			args[9] = vx[i];
			args[10] = vy[i];
			args[11] = vz[i];
			for (Int_t j = 0; j < 14; j++)
			{			
				args[12+j*3+0] = hitX[j][i];
				args[12+j*3+1] = hitY[j][i];
				args[12+j*3+2] = hitZ[j][i];
			}
			args[12+14*3+0] = dHLTsign[i];
			args[12+14*3+1] = dHLTpx[i];
			args[12+14*3+2] = dHLTpy[i];
			args[12+14*3+3] = dHLTpz[i];
			args[12+14*3+4] = dHLTfit[i];
			for (Int_t j = 0; j < 14; j++)
			{			
				args[12+14*3+5+j*3+0] = dHLThitX[j][i];
				args[12+14*3+5+j*3+1] = dHLThitY[j][i];
				args[12+14*3+5+j*3+2] = dHLThitZ[j][i];
			}
			table->Fill(args);
		}
		
	} // event loop

	TFile file("tracktable.root", "RECREATE");
	file.cd();
	table->Write(table->GetName(), TObject::kOverwrite);

	cout << "Done." << endl;
}
Commit	Line	Data
aa14645f	1	/**************************************************************************
	2	* This file is property of and copyright by the ALICE HLT Project *
	3	* All rights reserved. *
	4	* *
	5	* Primary Authors: *
	6	* Artur Szostak <artursz@iafrica.com> *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	// Simple macro to generate N-tuple for performance analysis.
	18	// This macro must be compiled and run like so from within the aliroot command
	19	// prompt:
b0635849	20	// root [0] gSystem->Load("libAliHLTMUON");
d8cf7391	21	// root [0] .x MakeHitsTable.C+
aa14645f	22
	23	#include "TVector3.h"
	24	#include "TSystem.h"
	25	#include "TROOT.h"
	26	#include "TParticle.h"
	27	#include "TArrayF.h"
	28	#include "TStopwatch.h"
	29	#include "TNtuple.h"
	30	#include "TFile.h"
	31	#include "TError.h"
af5746f2	32	#include "TMath.h"
aa14645f	33
846e6f41	34	#include "AliCDBManager.h"
aa14645f	35	#include "AliLoader.h"
	36	#include "AliRunLoader.h"
	37	#include "AliStack.h"
	38	#include "AliRun.h"
	39	#include "AliMUON.h"
	40	#include "AliMUONMCDataInterface.h"
	41	#include "AliMUONHit.h"
3b376c57	42	#include "AliHLTMUONEvent.h"
aa14645f	43	#include "AliHLTMUONRootifierComponent.h"
	44	#include "AliHLTMUONMansoTrack.h"
	45
	46	#include <cstdlib>
	47	#include <vector>
	48	#include <iostream>
	49	using std::cout;
	50	using std::cerr;
	51	using std::endl;
	52
	53
	54	void MakeTrackTable(
	55	Int_t firstEvent = 0,
	56	Int_t lastEvent = -1,
d8cf7391	57	const char* dHLToutputfile = "output.root",
af5746f2	58	Float_t maxSigma = 4., // 4 standard deviations
	59	Float_t sigmaX = 0.1, // 1 mm resolution
	60	Float_t sigmaY = 0.01, // 100 micron resolution
	61	Float_t sigmaZ = 0.02, // 200 microns resolution
	62	Float_t sigmaXtrg = 0.5, // 5 mm resolution
	63	Float_t sigmaYtrg = 0.5, // 5 mm resolution
	64	Float_t sigmaZtrg = 0.02 // 2 microns resolution
aa14645f	65	)
aa14645f	66	{
846e6f41	67	// Setup the CDB default storage and run number if nothing was set.
	68	AliCDBManager* cdbManager = AliCDBManager::Instance();
	69	if (cdbManager == NULL)
	70	{
	71	cerr << "ERROR: Global CDB manager object does not exist." << endl;
	72	return;
	73	}
	74	if (cdbManager->GetDefaultStorage() == NULL)
	75	{
162637e4	76	cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
846e6f41	77	}
	78	if (cdbManager->GetRun() == -1)
	79	{
	80	cdbManager->SetRun(0);
	81	}
	82
4070f709	83	gSystem->Load("libAliHLTMUON");
aa14645f	84	// Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.
	85
	86	TString fieldnames = "event:isprimary:hastrack:cantrigger:pdgcode:sign:px:py:pz";
	87	for (Int_t i = 0; i <= 14; i++)
	88	{
	89	fieldnames += ":x";
	90	fieldnames += i;
	91	fieldnames += ":y";
	92	fieldnames += i;
	93	fieldnames += ":z";
	94	fieldnames += i;
	95	}
	96	fieldnames += ":dHLTsign:dHLTpx:dHLTpy:dHLTpz:dHLTfit";
	97	for (Int_t i = 1; i <= 14; i++)
	98	{
	99	fieldnames += ":dHLTx";
	100	fieldnames += i;
	101	fieldnames += ":dHLTy";
	102	fieldnames += i;
	103	fieldnames += ":dHLTz";
	104	fieldnames += i;
	105	}
	106	TNtuple* table = new TNtuple("tracktable", "Table of tracks.", fieldnames);
	107
	108	TFile dHLTfile(dHLToutputfile, "READ");
	109
	110	AliMUONMCDataInterface data;
	111	AliStack* stack;
	112
	113	Int_t nEvents = data.NumberOfEvents();
	114	if (lastEvent < 0) lastEvent = nEvents-1;
	115	for (Int_t event = firstEvent; event <= lastEvent; event++)
	116	{
	117	cout << "Processing event: " << event;
	118
	119	char buf[1024];
	120	char* str = &buf[0];
	121	sprintf(str, "AliHLTMUONEvent;%d", event+1);
	122	AliHLTMUONEvent* dHLTevent = dynamic_cast<AliHLTMUONEvent*>( dHLTfile.Get(str) );
	123	if (dHLTevent == NULL)
	124	{
	125	cout << endl;
	126	cerr << "ERROR: Could not find " << str << " in " << dHLToutputfile << endl;
	127	continue;
	128	}
	129
	130	data.GetEvent(event);
	131
	132	stack = data.Stack(event);
	133	if (stack == NULL)
	134	{
	135	cout << endl;
	136	cerr << "ERROR: Stack is NULL" << endl;
	137	continue;
	138	}
	139
	140	Int_t trackcount = data.NumberOfParticles();
	141	cout << "\ttrackcount = " << trackcount << endl;
	142
	143	TArrayF isPrimary(trackcount);
	144	TArrayF hasTrack(trackcount);
	145	TArrayI hitsInTrigger(trackcount);
	146	TArrayF pdgcode(trackcount);
	147	TArrayF sign(trackcount);
148	TArrayF px(trackcount);
149	TArrayF py(trackcount);
150	TArrayF pz(trackcount);
151	TArrayF vx(trackcount);
152	TArrayF vy(trackcount);
153	TArrayF vz(trackcount);
154	TArrayF dHLTsign(trackcount);
155	TArrayF dHLTpx(trackcount);
156	TArrayF dHLTpy(trackcount);
157	TArrayF dHLTpz(trackcount);
158	TArrayF dHLTfit(trackcount);
159
160	//cout << "Fill particle data" << endl;
161	for (Int_t i = 0; i < trackcount; i++)
162	{
163	TParticle* p = data.Particle(i);
164	//isPrimary[i] = (i < stack->GetNprimary()) ? 1 : 0;
165	isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
166	hasTrack[i] = 0;
167	hitsInTrigger[i] = 0;
168	pdgcode[i] = p->GetPdgCode();
169	if (p->GetPDG() != NULL)
170	{
171	Double_t charge = p->GetPDG()->Charge();
172	if (charge > 0)
173	sign[i] = 1;
174	else if (charge < 0)
175	sign[i] = -1;
176	else
177	sign[i] = 0;
178	}
179	else
180	{
181	sign[i] = 0;
182	cerr << "ERROR: Could not get the PDG information." << endl;
183	}
184	px[i] = p->Px();
185	py[i] = p->Py();
186	pz[i] = p->Pz();
187	vx[i] = p->Vx();
188	vy[i] = p->Vy();
189	vz[i] = p->Vz();
190	dHLTsign[i] = 0;
191	dHLTpx[i] = 0;
192	dHLTpy[i] = 0;
193	dHLTpz[i] = 0;
194	dHLTfit[i] = -1;
195	}
196
197	TArrayF hitX[14];
198	TArrayF hitY[14];
199	TArrayF hitZ[14];
200	TArrayF dHLThitX[14];
201	TArrayF dHLThitY[14];
202	TArrayF dHLThitZ[14];
203	for (Int_t i = 0; i < 14; i++)
204	{
205	hitX[i].Set(trackcount);
206	hitY[i].Set(trackcount);
207	hitZ[i].Set(trackcount);
208	dHLThitX[i].Set(trackcount);
209	dHLThitY[i].Set(trackcount);
210	dHLThitZ[i].Set(trackcount);
211	for (Int_t j = 0; j < trackcount; j++)
212	{
213	hitX[i][j] = 0;
214	hitY[i][j] = 0;
215	hitZ[i][j] = 0;
216	dHLThitX[i][j] = 0;
217	dHLThitY[i][j] = 0;
218	dHLThitZ[i][j] = 0;
219	}
220	}
221
222	//cout << "Fill hits" << endl;
223	for (Int_t i = 0; i < data.NumberOfTracks(); i++)
224	for (Int_t j = 0; j < data.NumberOfHits(i); j++)
225	{
226	AliMUONHit* hit = data.Hit(i, j);
227	Int_t chamber = hit->Chamber() - 1;
228	if (chamber < 0 \|\| chamber >= 14)
229	{
230	cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
231	continue;
232	}
233
234	//cout << "hit->Track() = " << hit->Track() << endl;
235	if (hit->Track() < 0 \|\| hit->Track() >= trackcount)
236	{
237	cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
238	<< trackcount << ")" << endl;
239	continue;
240	}
241	Int_t particleindex = hit->Track();
242
243	hitX[chamber][particleindex] = hit->X();
244	hitY[chamber][particleindex] = hit->Y();
245	hitZ[chamber][particleindex] = hit->Z();
246
247	if (chamber == 0)
248	{
249	hasTrack[particleindex] = 1;
250	}
251	if (chamber >= 10)
252	{
253	hitsInTrigger[particleindex]++;
254	}
255	} // hits loop
256
257	// Ok now go through the dHLT information for this event and correlate
258	// the dHLT tracks to the kine data.
259	// This is done as follows: take a dHLT track. Then for each kine track calculate
260	// the distance between each hit of the kine track and corresponding dHLT track.
261	// Then sum up the distances and divide by the number of hits matched.
262	// This gives us a fit quality /Chi2 type parameter. The kine track that
263	// has the smallest sum of differences is then chosen as the kine track
264	// correlating to the current dHLT track being processed.
265	for (Int_t i = 0; i < dHLTevent->Array().GetEntriesFast(); i++)
266	{
267	if (dHLTevent->Array().At(i)->IsA() != AliHLTMUONMansoTrack::Class())
268	continue;
269	AliHLTMUONMansoTrack* mtrack = static_cast<AliHLTMUONMansoTrack*>(dHLTevent->Array().At(i));
270
271	// Find the best fitting kine track.
272	Double_t bestFitQuality = 1e30;
273	Int_t bestFitTrack = -1;
274	for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++)
275	{
276	Double_t sumOfDiffs = 0;
277	Double_t hitsMatched = 0;
278	for (Int_t j = 7; j <= 10; j++)
279	{
280	const AliHLTMUONRecHit* hit = mtrack->Hit(j);
281	if (hit == NULL) continue;
282	if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
283	continue;
284
285	TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
286	TVector3 diff = hV - hit->Coordinate();
af5746f2	287	Double_t diffX = diff.X() / sigmaX;
	288	Double_t diffY = diff.Y() / sigmaY;
	289	Double_t diffZ = diff.Z() / sigmaZ;
	290	if (diffX > maxSigma) continue;
	291	if (diffY > maxSigma) continue;
	292	if (diffZ > maxSigma) continue;
	293	sumOfDiffs += diffXdiffX + diffYdiffY + diffZ*diffZ;
aa14645f	294	hitsMatched++;
	295	}
	296	if (mtrack->TriggerRecord() != NULL)
	297	{
	298	for (Int_t j = 11; j <= 14; j++)
	299	{
	300	const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
	301	if (hit == TVector3(0,0,0)) continue;
	302	if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
	303	continue;
	304
	305	TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
	306	TVector3 diff = hV - hit;
af5746f2	307	Double_t diffX = diff.X() / sigmaXtrg;
	308	Double_t diffY = diff.Y() / sigmaYtrg;
	309	Double_t diffZ = diff.Z() / sigmaZtrg;
	310	if (diffX > maxSigma) continue;
	311	if (diffY > maxSigma) continue;
	312	if (diffZ > maxSigma) continue;
	313	sumOfDiffs += diffXdiffX + diffYdiffY + diffZ*diffZ;
aa14645f	314	hitsMatched++;
	315	}
	316	}
	317
	318	// Now check the fit quality.
	319	if (hitsMatched <= 0) continue;
af5746f2	320	Double_t fitQuality = TMath::Sqrt(sumOfDiffs) / hitsMatched;
aa14645f	321	if (fitQuality < bestFitQuality)
	322	{
	323	bestFitQuality = fitQuality;
	324	bestFitTrack = kinetrack;
	325	}
	326	}
	327
	328	if (bestFitTrack < 0)
	329	{
	330	// No track was matched so we need to add a fake track record.
	331	Float_t args[12+143+5+143];
	332	args[0] = event;
	333	args[1] = -1;
	334	args[2] = 0;
	335	args[3] = 0;
	336	args[4] = 0;
	337	args[5] = 0;
	338	args[6] = 0;
	339	args[7] = 0;
	340	args[8] = 0;
	341	args[9] = 0;
	342	args[10] = 0;
	343	args[11] = 0;
	344	for (Int_t j = 0; j < 14; j++)
	345	{
	346	args[12+j*3+0] = 0;
	347	args[12+j*3+1] = 0;
	348	args[12+j*3+2] = 0;
	349	}
	350	args[12+14*3+0] = mtrack->Sign();
	351	args[12+14*3+1] = mtrack->Px();
	352	args[12+14*3+2] = mtrack->Py();
	353	args[12+14*3+3] = mtrack->Pz();
	354	args[12+14*3+4] = -1;
	355	for (Int_t j = 0; j < 14; j++)
	356	{
	357	args[12+143+5+j3+0] = 0;
	358	args[12+143+5+j3+1] = 0;
	359	args[12+143+5+j3+2] = 0;
	360	}
	361	for (Int_t j = 6; j < 10; j++)
	362	{
	363	const AliHLTMUONRecHit* hit = mtrack->Hit(j+1);
	364	if (hit == NULL) continue;
	365	args[12+143+5+j3+0] = hit->X();
	366	args[12+143+5+j3+1] = hit->Y();
	367	args[12+143+5+j3+2] = hit->Z();
	368	}
	369	if (mtrack->TriggerRecord() != NULL)
	370	{
	371	for (Int_t j = 10; j < 14; j++)
	372	{
	373	const TVector3& hit = mtrack->TriggerRecord()->Hit(j+1);
	374	args[12+143+5+j3+0] = hit.X();
	375	args[12+143+5+j3+1] = hit.Y();
	376	args[12+143+5+j3+2] = hit.Z();
	377	}
	378	}
	379	table->Fill(args);
	380	}
	381	else
	382	{
	383	// Fill the details about the dHLT track to the best fitting track info.
	384	dHLTsign[bestFitTrack] = mtrack->Sign();
385	dHLTpx[bestFitTrack] = mtrack->Px();
386	dHLTpy[bestFitTrack] = mtrack->Py();
387	dHLTpz[bestFitTrack] = mtrack->Pz();
388	dHLTfit[bestFitTrack] = bestFitQuality;
389	for (Int_t j = 7; j <= 10; j++)
390	{
391	const AliHLTMUONRecHit* hit = mtrack->Hit(j);
392	if (hit == NULL) continue;
393	dHLThitX[j-1][bestFitTrack] = hit->X();
394	dHLThitY[j-1][bestFitTrack] = hit->Y();
395	dHLThitZ[j-1][bestFitTrack] = hit->Z();
396	}
397	if (mtrack->TriggerRecord() != NULL)
398	{
399	for (Int_t j = 11; j <= 14; j++)
400	{
401	const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
402	if (hit == TVector3(0,0,0)) continue;
403	dHLThitX[j-1][bestFitTrack] = hit.X();
404	dHLThitY[j-1][bestFitTrack] = hit.Y();
405	dHLThitZ[j-1][bestFitTrack] = hit.Z();
406	}
407	}
408	}
409	}
410
411	//cout << "Fill table" << endl;
412	for (Int_t i = 0; i < trackcount; i++)
413	{
414	Float_t args[12+143+5+143];
415	args[0] = event;
416	args[1] = isPrimary[i];
417	args[2] = hasTrack[i];
418	args[3] = (hitsInTrigger[i] > 2) ? 1 : 0;
419	args[4] = pdgcode[i];
420	args[5] = sign[i];
421	args[6] = px[i];
422	args[7] = py[i];
423	args[8] = pz[i];
424	args[9] = vx[i];
425	args[10] = vy[i];
426	args[11] = vz[i];
427	for (Int_t j = 0; j < 14; j++)
428	{
429	args[12+j*3+0] = hitX[j][i];
430	args[12+j*3+1] = hitY[j][i];
431	args[12+j*3+2] = hitZ[j][i];
432	}
433	args[12+14*3+0] = dHLTsign[i];
434	args[12+14*3+1] = dHLTpx[i];
435	args[12+14*3+2] = dHLTpy[i];
436	args[12+14*3+3] = dHLTpz[i];
437	args[12+14*3+4] = dHLTfit[i];
438	for (Int_t j = 0; j < 14; j++)
439	{
440	args[12+143+5+j3+0] = dHLThitX[j][i];
441	args[12+143+5+j3+1] = dHLThitY[j][i];
442	args[12+143+5+j3+2] = dHLThitZ[j][i];
443	}
444	table->Fill(args);
445	}
446
447	} // event loop
448
449	TFile file("tracktable.root", "RECREATE");
450	file.cd();
451	table->Write(table->GetName(), TObject::kOverwrite);
452
453	cout << "Done." << endl;
454	}