[u/mrichter/AliRoot.git] / ACORDE / macros / AnalysisMacros / Local / AliAnalysisTaskAcorde.cxx

/////////////////////////////////////////////////////////////////////////////
//
// AliAnalysisTaskAcorde class
//
// Description:
//
//	Reads the information of ACORDE-ESD
//	Also it is implemented a simple matching between tracks
//	to look for the extrapolation of them to ACORDE modules
//
//	Create some histograms and one tree with the information
//	of the matching tracks
//
//  Author: Mario Rodríguez Cahuantzi
//		<mario.rocah@gmail.com>
//		<mrodrigu@mail.cern.ch>
//
//  Created: June 30th. 2010 @ FCFM - BUAP, Puebla, MX
//  Last update: created
//
/////////////////////////////////////////////////////////////////////////////
#include "TChain.h"
#include "TTree.h"
#include "TH2F.h"
#include "TH1F.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"

#include "AliAnalysisTaskAcorde.h"

#include "TMath.h"
#include "TArrayI.h"
#include "TDatabasePDG.h"
#include "TVectorD.h"

#include "AliESDtrack.h"
#include "AliTracker.h"
#include "TFile.h"
#include "TVectorD.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
ClassImp(AliAnalysisTaskAcorde)
//________________________________________________________________________
AliAnalysisTaskAcorde::AliAnalysisTaskAcorde(const char *name) 
  : AliAnalysisTask(name, ""), 
	fESD(0),
	cosmicTree(0),
	nTracks(0),
	nMatchTracks(0),
	minTPCclusters(30),
	minTrackDist(1000.),
	minCutDir(0.95),
  	xAco(0),
  	yAco(0),
  	zAco(0),
	trigger(0),
	ActiveTriggerDetector(0),
	nSingleTracks(0),
	nMatchedTracks(0),
	histo(0),
 	ntracks(0),
	acordeHitsAll(0),
	acordeMultiAll(0),
	acordeHitsTPC(0),
	acordeMultiTPC(0),
	nTracksMatched(0),
	fTracksToAcorde(0)

{
	// Constructor of class

  	DefineInput(0, TChain::Class());
    	DefineOutput(0,TTree::Class());
	DefineOutput(1,TList::Class());
}

//________________________________________________________________________
AliAnalysisTaskAcorde::~AliAnalysisTaskAcorde()
{
	// destructor --> avoid watchdog mails?
	
	delete fESD;
	delete histo;
	delete cosmicTree;
}
//________________________________________________________________________
void AliAnalysisTaskAcorde::ConnectInputData(Option_t *) 
{
	// Connect ESD or AOD here
	// Called once

	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  	if (!tree) 
	{
    		Printf("ERROR: Could not read chain from input slot 0");
  	}else 
	{

    		AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
    		if (!esdH) 
		{
      			Printf("ERROR: Could not get ESDInputHandler");
    		}else fESD = esdH->GetEvent();
  	}
}

//________________________________________________________________________
void AliAnalysisTaskAcorde::CreateOutputObjects()
{

	// creates one TList with some histograms

	histo = new TList();

	nTracksMatched = new TH1F("nTracksMatched","nTracksMatched from matching algorithm implemented",300,1,300);
	ntracks = new TH1F("ntracks","nTracks distribution",500,1,500);
	acordeHitsAll = new TH1F("acordeHitsAll","Hits of ACORDE from ESD",60,-0.5,59.5);
	acordeHitsTPC = new TH1F("acordeHitsAllTPC","Hits of ACORDE from ESD (if ntracks>0)",61,0,60);
	acordeMultiAll = new TH1F("acordeMultiAll","Multiplicity of ACORDE modules from ESD",60,-0.5,59.5);
	acordeMultiTPC = new TH1F("acordeMultiAllTPC","Multiplicity of ACORDE modules from ESD (id ntracks>0)",61,0,60);
	fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);

	histo->Add(nTracksMatched);
	histo->Add(ntracks);
	histo->Add(acordeHitsAll);
	histo->Add(acordeHitsTPC);
	histo->Add(acordeMultiAll);
	histo->Add(acordeMultiTPC);
	histo->Add(fTracksToAcorde);

	// Create Tree branches
  	// Called just once

  	cosmicTree = new TTree("cosmicTree","cosmicTreeMRC");
	cosmicTree->Branch("nTracks",&nTracks,"nTracks/I");
	cosmicTree->Branch("nMatchTracks",&nMatchTracks,"nMatchTracks/I");
	cosmicTree->Branch("xAco",&xAco,"xAco/F");
	cosmicTree->Branch("yAco",&yAco,"yAco/F");
	cosmicTree->Branch("zAco",&zAco,"zAco/F");
	cosmicTree->Branch("trigger",&trigger,"trigger/I");

}


void AliAnalysisTaskAcorde::Exec(Option_t *) 
{
	// Main loop
  	// Called for each event

	Int_t counterOfMuons = 0;


  	if (!fESD) 
	{
    		Printf("ERROR: fESD not available");
    		return;
  	}


	// Pointer to the information of ACORDE detector

	AliESDACORDE *acordeESD = fESD->GetACORDEData();

	Int_t contMulti = 0;
	Int_t contMultiTPC = 0;

	for(Int_t imod=0;imod<60;imod++)
	{
		if (acordeESD->GetHitChannel(imod)) 
		{
			acordeHitsAll->Fill(imod);
			contMulti++;
		}

	}acordeMultiAll->Fill(contMulti);

	for(Int_t imod=0;imod<60;imod++)
	{
		if (acordeESD->GetHitChannel(imod))
		{
			acordeHitsTPC->Fill(imod);
			contMultiTPC++;
		}
	}acordeMultiTPC->Fill(contMultiTPC);


	// Assingment of the type of Trigger to cosmicTree

	ActiveTriggerDetector = fESD->GetFiredTriggerClasses();

	if (ActiveTriggerDetector.Contains("C0SCO")) trigger=0;
	if (ActiveTriggerDetector.Contains("C0SH2")) trigger=1;
  	if (ActiveTriggerDetector.Contains("C0AMU")) trigger=2;
  	if (ActiveTriggerDetector.Contains("C0LSR")) trigger=3;
  	if (ActiveTriggerDetector.Contains("C0SCO1")) trigger=4;
  	if (ActiveTriggerDetector.Contains("C0OBE")) trigger=5;
  	if (ActiveTriggerDetector.Contains("C0PCO")) trigger=6;
  	if (ActiveTriggerDetector.Contains("C0OB3")) trigger=7;
  	if (ActiveTriggerDetector.Contains("C0OCP")) trigger=8;
    	if (ActiveTriggerDetector.Contains("C0ASL")) trigger=9;
	

	// Begins the matching analysis between tracks

	nTracks = fESD->GetNumberOfTracks();
	if (nTracks<=0) return;
	ntracks->Fill(nTracks);

	fPair = new TArrayI(nTracks);
	nSingleTracks=0;
	nMatchedTracks=0;
	Int_t muons = 0;
	TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
	for (Int_t iTrack=0; iTrack<nTracks; iTrack++)
	{

		// Special case: nTracks == 1
		if (nTracks ==1)
		{
			AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
			if (!singleTrack || !singleTrack->GetOuterParam() || !singleTrack->GetInnerParam()) continue;
			if (singleTrack->GetTPCNcls()<minTPCclusters) continue;
			Float_t vertX = singleTrack->Xv();
			Float_t vertZ = singleTrack->Zv();
			if (TMath::Abs(vertX)>165 || TMath::Abs(vertZ)>245) continue; // far away from the VERTEX
			if (singleTrack->E()< 1) continue; // less than 1 GeV 
			if (singleTrack->P()< 1) continue; // less than 1 GeV 
			nMatchTracks=0;
			nSingleTracks = nTracks;
			muons = 1;

		} // nTracks ==1
		else if (nTracks>1) // nTracks > 1
		{

			AliESDtrack *track0 = fESD->GetTrack(iTrack);
			(*fPair)[iTrack]=-1;
			if (!track0 || !track0->GetOuterParam() || !track0->GetInnerParam()) continue;

			Double_t track0_dir[3];
			track0->GetDirection(track0_dir);
			Float_t minDist = minTrackDist;	

			const AliExternalTrackParam * trackIn = track0->GetInnerParam();
			const AliExternalTrackParam * trackOut = track0->GetOuterParam();

			if (!trackIn || !trackOut) continue;
			if (nTracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // Filter TPC-Laser

			if (track0->GetTPCNcls()<minTPCclusters) continue;

			for (Int_t jTrack=iTrack; jTrack<nTracks; ++jTrack)
			{
				AliESDtrack *track1 = fESD->GetTrack(jTrack);
				if (!track1 || !track1->GetOuterParam() || !track1->GetInnerParam()) continue;
				Double_t track1_dir[3];
				track1->GetDirection(track1_dir);
				if (track1->GetTPCNcls()<minTPCclusters) continue;
				Float_t direction = track0_dir[0]*track1_dir[0] + track0_dir[1]*track1_dir[1] + track0_dir[2]*track1_dir[2];
				if (TMath::Abs(direction)<minCutDir) continue; // direction vector product
				Float_t dvertex0[3];
				Float_t dvertex1[3];

				dvertex0[0] = track0->Xv();
				dvertex0[1] = track0->Yv();
				dvertex0[2] = track0->Zv();

				dvertex1[0] = track1->Xv();
				dvertex1[1] = track1->Yv();
				dvertex1[2] = track1->Zv();

				if (TMath::Abs(dvertex0[0])>165 || TMath::Abs(dvertex0[2])>245 || TMath::Abs(dvertex1[0])>165 || TMath::Abs(dvertex0[2])>245) continue;
				// calculate the distance between track0 and track1
				
				Float_t dy = track0->GetY()+track1->GetY();
				Float_t sy2 = track0->GetSigmaY2()+track1->GetSigmaY2();
				Float_t dphi = track0->GetAlpha()-track1->GetAlpha()-TMath::Pi();
				Float_t sphi2 = track0->GetSigmaSnp2()+track1->GetSigmaSnp2();
				Float_t dtheta = track0->GetTgl()+track1->GetTgl();
				Float_t stheta2 = track0->GetSigmaTgl2()+track1->GetSigmaTgl2();
				Float_t normDist = TMath::Sqrt(dy*dy/sy2+dphi+dphi/sphi2+dtheta*dtheta/stheta2);
				if (normDist>minDist) continue;

				Float_t distTracks = TMath::Sqrt((dvertex0[0]-dvertex1[0])*(dvertex0[0]-dvertex1[0]) + (dvertex0[1]-dvertex1[1])*(dvertex0[1]-dvertex1[1]) + (dvertex0[2]-dvertex1[2])*(dvertex0[2]-dvertex1[2]) );
				if (distTracks>2.5) continue;


				minDist = normDist;

				// after all the cuts we should have only tracks that can be matched

				(*fPair)[iTrack] = jTrack;


			}	
		} // nTracks > 1
	

	} // Loop over all the tracks


	// Matching tracks
	
	for (Int_t itrack=0;itrack<nTracks;itrack++)
	{
		AliESDtrack *UpTrack = 0;
		AliESDtrack *DownTrack = 0;

		AliESDtrack *Track0 = fESD->GetTrack(itrack);
		if (!Track0 || !Track0->GetOuterParam() || !Track0->GetInnerParam() ) continue;

		Double_t mxyz[3];
		Track0->GetOuterParam()->GetXYZ(mxyz);
		
		UpTrack = Track0;
	

		if ((*fPair)[itrack]>0)
		{
			AliESDtrack *Track1 = fESD->GetTrack((*fPair)[itrack]);
			if (!Track1 || !Track1->GetOuterParam() || !Track1->GetInnerParam()) continue;
			Float_t dvertex0[3];
			Float_t dvertex1[3];

			dvertex0[0] = Track0->Xv();
			dvertex0[1] = Track0->Yv();
			dvertex0[2] = Track0->Zv();

			dvertex1[0] = Track1->Xv();
			dvertex1[1] = Track1->Yv();
			dvertex1[2] = Track1->Zv();

			Float_t distTracks = TMath::Sqrt((dvertex0[0]-dvertex1[0])*(dvertex0[0]-dvertex1[0]) + (dvertex0[1]-dvertex1[1])*(dvertex0[1]-dvertex1[1]) + (dvertex0[2]-dvertex1[2])*(dvertex0[2]-dvertex1[2]) );
			if (distTracks>2.5) continue;

			if (Track1->GetOuterParam())
			{
				Double_t nxyz[3];
				Track1->GetOuterParam()->GetXYZ(nxyz);
				if (nxyz[1]>mxyz[1])
				{
					UpTrack = Track1;
				}else DownTrack = Track1;
			}


			// Here we propagate the Up-Tracks to ACORDE


			const Double_t kRL3 = 510; // radious of L3 magnet
			const Double_t kxAcorde = 850.; // radios of "ACORDE" above the magnet

			Double_t agxyz[3];
			

			AliExternalTrackParam *upper = (AliExternalTrackParam *)(UpTrack->GetOuterParam()->Clone());

			Bool_t isOk = upper->PropagateTo(kRL3,fESD->GetMagneticField());
			upper->GetXYZ(agxyz);
			if (agxyz[1]<0) continue;

			upper->GetXYZ(agxyz);
			Double_t galpha = TMath::ATan2(agxyz[1],agxyz[0]);
			Double_t alpha = galpha;
			galpha*=180/TMath::Pi();

			if (galpha<45.) alpha = TMath::Pi()/8;
			if (galpha>135.) alpha = TMath::Pi()*(7/8);
			if (galpha>45 && galpha<135.) alpha = TMath::Pi()/2;

			if (isOk) upper->Rotate(alpha);

			if (isOk) isOk = upper->PropagateTo(kxAcorde,0);

			TVectorD rgxyz(3);
			AliExternalTrackParam *param = (AliExternalTrackParam *)upper;
			param->GetXYZ(rgxyz.GetMatrixArray());
			xAco = rgxyz[0];
			yAco = rgxyz[1];
			zAco = rgxyz[2];
			fTracksToAcorde->Fill(xAco,zAco);
			// Count how many tracks have been matched
			nMatchedTracks++;
		} else nSingleTracks++;
	}

	if ( (nMatchedTracks+nSingleTracks) <= nTracks ) muons=nMatchedTracks+nSingleTracks;
	nTracksMatched->Fill(muons);

	nMatchTracks = nSingleTracks + nMatchedTracks;
	
	// Fills the tree

	cosmicTree->Fill();
	
	// Post output data.

	PostData(0,cosmicTree);
	PostData(1,histo);
}      

//________________________________________________________________________
void AliAnalysisTaskAcorde::Terminate(Option_t *) 
{

//	AliDebug(1,"Do nothig in Terminate");

	//nPart->Draw();
        TCanvas *c1 = new TCanvas("c1","TPC-resolution for P reconstruction .... MRC");

        c1->Divide(2,2);
	c1->cd(1);
	acordeHitsAll->Draw();
	c1->cd(2)->SetLogy();
	acordeMultiAll->Draw();
	c1->cd(3)->SetLogy();
	//nTracksMatched->Draw();
	ntracks->Draw();
	c1->cd(4);
	fTracksToAcorde->Draw();
	
	
}
Commit	Line	Data
b5f2577c	1	/////////////////////////////////////////////////////////////////////////////
	2	//
	3	// AliAnalysisTaskAcorde class
	4	//
	5	// Description:
	6	//
	7	// Reads the information of ACORDE-ESD
	8	// Also it is implemented a simple matching between tracks
	9	// to look for the extrapolation of them to ACORDE modules
	10	//
	11	// Create some histograms and one tree with the information
	12	// of the matching tracks
	13	//
	14	// Author: Mario Rodríguez Cahuantzi
	15	// <mario.rocah@gmail.com>
	16	// <mrodrigu@mail.cern.ch>
	17	//
	18	// Created: June 30th. 2010 @ FCFM - BUAP, Puebla, MX
	19	// Last update: created
	20	//
	21	/////////////////////////////////////////////////////////////////////////////
	22	#include "TChain.h"
	23	#include "TTree.h"
	24	#include "TH2F.h"
	25	#include "TH1F.h"
	26	#include "AliAnalysisTask.h"
	27	#include "AliAnalysisManager.h"
	28
	29	#include "AliESDEvent.h"
	30	#include "AliESDInputHandler.h"
	31
	32	#include "AliAnalysisTaskAcorde.h"
	33
	34	#include "TMath.h"
	35	#include "TArrayI.h"
	36	#include "TDatabasePDG.h"
	37	#include "TVectorD.h"
	38
	39	#include "AliESDtrack.h"
	40	#include "AliTracker.h"
	41	#include "TFile.h"
	42	#include "TVectorD.h"
	43	#include "TStyle.h"
	44	#include "TCanvas.h"
	45	#include "TLegend.h"
	46	ClassImp(AliAnalysisTaskAcorde)
	47	//________________________________________________________________________
	48	AliAnalysisTaskAcorde::AliAnalysisTaskAcorde(const char *name)
	49	: AliAnalysisTask(name, ""),
	50	fESD(0),
	51	cosmicTree(0),
	52	nTracks(0),
	53	nMatchTracks(0),
	54	minTPCclusters(30),
	55	minTrackDist(1000.),
	56	minCutDir(0.95),
	57	xAco(0),
	58	yAco(0),
	59	zAco(0),
	60	trigger(0),
	61	ActiveTriggerDetector(0),
	62	nSingleTracks(0),
	63	nMatchedTracks(0),
	64	histo(0),
65	ntracks(0),
66	acordeHitsAll(0),
67	acordeMultiAll(0),
68	acordeHitsTPC(0),
69	acordeMultiTPC(0),
70	nTracksMatched(0),
71	fTracksToAcorde(0)
72
73	{
74	// Constructor of class
75
76	DefineInput(0, TChain::Class());
77	DefineOutput(0,TTree::Class());
78	DefineOutput(1,TList::Class());
79	}
80
81	//________________________________________________________________________
82	AliAnalysisTaskAcorde::~AliAnalysisTaskAcorde()
83	{
84	// destructor --> avoid watchdog mails?
85
86	delete fESD;
87	delete histo;
88	delete cosmicTree;
89	}
90	//________________________________________________________________________
91	void AliAnalysisTaskAcorde::ConnectInputData(Option_t *)
92	{
93	// Connect ESD or AOD here
94	// Called once
95
96	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
97	if (!tree)
98	{
99	Printf("ERROR: Could not read chain from input slot 0");
100	}else
101	{
102
103	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
104	if (!esdH)
105	{
106	Printf("ERROR: Could not get ESDInputHandler");
107	}else fESD = esdH->GetEvent();
108	}
109	}
110
111	//________________________________________________________________________
112	void AliAnalysisTaskAcorde::CreateOutputObjects()
113	{
114
115	// creates one TList with some histograms
116
117	histo = new TList();
118
119	nTracksMatched = new TH1F("nTracksMatched","nTracksMatched from matching algorithm implemented",300,1,300);
120	ntracks = new TH1F("ntracks","nTracks distribution",500,1,500);
121	acordeHitsAll = new TH1F("acordeHitsAll","Hits of ACORDE from ESD",60,-0.5,59.5);
122	acordeHitsTPC = new TH1F("acordeHitsAllTPC","Hits of ACORDE from ESD (if ntracks>0)",61,0,60);
123	acordeMultiAll = new TH1F("acordeMultiAll","Multiplicity of ACORDE modules from ESD",60,-0.5,59.5);
124	acordeMultiTPC = new TH1F("acordeMultiAllTPC","Multiplicity of ACORDE modules from ESD (id ntracks>0)",61,0,60);
125	fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);
126
127	histo->Add(nTracksMatched);
128	histo->Add(ntracks);
129	histo->Add(acordeHitsAll);
130	histo->Add(acordeHitsTPC);
131	histo->Add(acordeMultiAll);
132	histo->Add(acordeMultiTPC);
133	histo->Add(fTracksToAcorde);
134
135	// Create Tree branches
136	// Called just once
137
138	cosmicTree = new TTree("cosmicTree","cosmicTreeMRC");
139	cosmicTree->Branch("nTracks",&nTracks,"nTracks/I");
140	cosmicTree->Branch("nMatchTracks",&nMatchTracks,"nMatchTracks/I");
141	cosmicTree->Branch("xAco",&xAco,"xAco/F");
142	cosmicTree->Branch("yAco",&yAco,"yAco/F");
143	cosmicTree->Branch("zAco",&zAco,"zAco/F");
144	cosmicTree->Branch("trigger",&trigger,"trigger/I");
145
146	}
147
148
149	void AliAnalysisTaskAcorde::Exec(Option_t *)
150	{
151	// Main loop
152	// Called for each event
153
154	Int_t counterOfMuons = 0;
155
156
157
158	if (!fESD)
159	{
160	Printf("ERROR: fESD not available");
161	return;
162	}
163
164
165
166	// Pointer to the information of ACORDE detector
167
168	AliESDACORDE *acordeESD = fESD->GetACORDEData();
169
170	Int_t contMulti = 0;
171	Int_t contMultiTPC = 0;
172
173	for(Int_t imod=0;imod<60;imod++)
174	{
175	if (acordeESD->GetHitChannel(imod))
176	{
177	acordeHitsAll->Fill(imod);
178	contMulti++;
179	}
180
181	}acordeMultiAll->Fill(contMulti);
182
183	for(Int_t imod=0;imod<60;imod++)
184	{
185	if (acordeESD->GetHitChannel(imod))
186	{
187	acordeHitsTPC->Fill(imod);
188	contMultiTPC++;
189	}
190	}acordeMultiTPC->Fill(contMultiTPC);
191
192
193	// Assingment of the type of Trigger to cosmicTree
194
195	ActiveTriggerDetector = fESD->GetFiredTriggerClasses();
196
197	if (ActiveTriggerDetector.Contains("C0SCO")) trigger=0;
198	if (ActiveTriggerDetector.Contains("C0SH2")) trigger=1;
199	if (ActiveTriggerDetector.Contains("C0AMU")) trigger=2;
200	if (ActiveTriggerDetector.Contains("C0LSR")) trigger=3;
201	if (ActiveTriggerDetector.Contains("C0SCO1")) trigger=4;
202	if (ActiveTriggerDetector.Contains("C0OBE")) trigger=5;
203	if (ActiveTriggerDetector.Contains("C0PCO")) trigger=6;
204	if (ActiveTriggerDetector.Contains("C0OB3")) trigger=7;
205	if (ActiveTriggerDetector.Contains("C0OCP")) trigger=8;
206	if (ActiveTriggerDetector.Contains("C0ASL")) trigger=9;
207
208
209	// Begins the matching analysis between tracks
210
211	nTracks = fESD->GetNumberOfTracks();
212	if (nTracks<=0) return;
213	ntracks->Fill(nTracks);
214
215	fPair = new TArrayI(nTracks);
216	nSingleTracks=0;
217	nMatchedTracks=0;
218	Int_t muons = 0;
219	TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
220	for (Int_t iTrack=0; iTrack<nTracks; iTrack++)
221	{
222
223	// Special case: nTracks == 1
224	if (nTracks ==1)
225	{
226	AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
227	if (!singleTrack \|\| !singleTrack->GetOuterParam() \|\| !singleTrack->GetInnerParam()) continue;
228	if (singleTrack->GetTPCNcls()<minTPCclusters) continue;
229	Float_t vertX = singleTrack->Xv();
230	Float_t vertZ = singleTrack->Zv();
231	if (TMath::Abs(vertX)>165 \|\| TMath::Abs(vertZ)>245) continue; // far away from the VERTEX
232	if (singleTrack->E()< 1) continue; // less than 1 GeV
233	if (singleTrack->P()< 1) continue; // less than 1 GeV
234	nMatchTracks=0;
235	nSingleTracks = nTracks;
236	muons = 1;
237
238	} // nTracks ==1
239	else if (nTracks>1) // nTracks > 1
240	{
241
242	AliESDtrack *track0 = fESD->GetTrack(iTrack);
243	(*fPair)[iTrack]=-1;
244	if (!track0 \|\| !track0->GetOuterParam() \|\| !track0->GetInnerParam()) continue;
245
246	Double_t track0_dir[3];
247	track0->GetDirection(track0_dir);
248	Float_t minDist = minTrackDist;
249
250	const AliExternalTrackParam * trackIn = track0->GetInnerParam();
251	const AliExternalTrackParam * trackOut = track0->GetOuterParam();
252
253	if (!trackIn \|\| !trackOut) continue;
254	if (nTracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // Filter TPC-Laser
255
256	if (track0->GetTPCNcls()<minTPCclusters) continue;
257
258	for (Int_t jTrack=iTrack; jTrack<nTracks; ++jTrack)
259	{
260	AliESDtrack *track1 = fESD->GetTrack(jTrack);
261	if (!track1 \|\| !track1->GetOuterParam() \|\| !track1->GetInnerParam()) continue;
262	Double_t track1_dir[3];
263	track1->GetDirection(track1_dir);
264	if (track1->GetTPCNcls()<minTPCclusters) continue;
265	Float_t direction = track0_dir[0]track1_dir[0] + track0_dir[1]track1_dir[1] + track0_dir[2]*track1_dir[2];
266	if (TMath::Abs(direction)<minCutDir) continue; // direction vector product
267	Float_t dvertex0[3];
268	Float_t dvertex1[3];
269
270	dvertex0[0] = track0->Xv();
271	dvertex0[1] = track0->Yv();
272	dvertex0[2] = track0->Zv();
273
274	dvertex1[0] = track1->Xv();
275	dvertex1[1] = track1->Yv();
276	dvertex1[2] = track1->Zv();
277
278	if (TMath::Abs(dvertex0[0])>165 \|\| TMath::Abs(dvertex0[2])>245 \|\| TMath::Abs(dvertex1[0])>165 \|\| TMath::Abs(dvertex0[2])>245) continue;
279	// calculate the distance between track0 and track1
280
281	Float_t dy = track0->GetY()+track1->GetY();
282	Float_t sy2 = track0->GetSigmaY2()+track1->GetSigmaY2();
283	Float_t dphi = track0->GetAlpha()-track1->GetAlpha()-TMath::Pi();
284	Float_t sphi2 = track0->GetSigmaSnp2()+track1->GetSigmaSnp2();
285	Float_t dtheta = track0->GetTgl()+track1->GetTgl();
286	Float_t stheta2 = track0->GetSigmaTgl2()+track1->GetSigmaTgl2();
287	Float_t normDist = TMath::Sqrt(dydy/sy2+dphi+dphi/sphi2+dthetadtheta/stheta2);
288	if (normDist>minDist) continue;
289
290	Float_t distTracks = TMath::Sqrt((dvertex0[0]-dvertex1[0])(dvertex0[0]-dvertex1[0]) + (dvertex0[1]-dvertex1[1])(dvertex0[1]-dvertex1[1]) + (dvertex0[2]-dvertex1[2])*(dvertex0[2]-dvertex1[2]) );
291	if (distTracks>2.5) continue;
292
293
294	minDist = normDist;
295
296	// after all the cuts we should have only tracks that can be matched
297
298	(*fPair)[iTrack] = jTrack;
299
300
301	}
302	} // nTracks > 1
303
304
305	} // Loop over all the tracks
306
307
308
309
310	// Matching tracks
311
312	for (Int_t itrack=0;itrack<nTracks;itrack++)
313	{
314	AliESDtrack *UpTrack = 0;
315	AliESDtrack *DownTrack = 0;
316
317	AliESDtrack *Track0 = fESD->GetTrack(itrack);
318	if (!Track0 \|\| !Track0->GetOuterParam() \|\| !Track0->GetInnerParam() ) continue;
319
320	Double_t mxyz[3];
321	Track0->GetOuterParam()->GetXYZ(mxyz);
322
323	UpTrack = Track0;
324
325
326	if ((*fPair)[itrack]>0)
327	{
328	AliESDtrack Track1 = fESD->GetTrack((fPair)[itrack]);
329	if (!Track1 \|\| !Track1->GetOuterParam() \|\| !Track1->GetInnerParam()) continue;
330	Float_t dvertex0[3];
331	Float_t dvertex1[3];
332
333	dvertex0[0] = Track0->Xv();
334	dvertex0[1] = Track0->Yv();
335	dvertex0[2] = Track0->Zv();
336
337	dvertex1[0] = Track1->Xv();
338	dvertex1[1] = Track1->Yv();
339	dvertex1[2] = Track1->Zv();
340
341	Float_t distTracks = TMath::Sqrt((dvertex0[0]-dvertex1[0])(dvertex0[0]-dvertex1[0]) + (dvertex0[1]-dvertex1[1])(dvertex0[1]-dvertex1[1]) + (dvertex0[2]-dvertex1[2])*(dvertex0[2]-dvertex1[2]) );
342	if (distTracks>2.5) continue;
343
344	if (Track1->GetOuterParam())
345	{
346	Double_t nxyz[3];
347	Track1->GetOuterParam()->GetXYZ(nxyz);
348	if (nxyz[1]>mxyz[1])
349	{
350	UpTrack = Track1;
351	}else DownTrack = Track1;
352	}
353
354
355
356	// Here we propagate the Up-Tracks to ACORDE
357
358
359	const Double_t kRL3 = 510; // radious of L3 magnet
360	const Double_t kxAcorde = 850.; // radios of "ACORDE" above the magnet
361
362	Double_t agxyz[3];
363
364
365	AliExternalTrackParam upper = (AliExternalTrackParam )(UpTrack->GetOuterParam()->Clone());
366
367	Bool_t isOk = upper->PropagateTo(kRL3,fESD->GetMagneticField());
368	upper->GetXYZ(agxyz);
369	if (agxyz[1]<0) continue;
370
371	upper->GetXYZ(agxyz);
372	Double_t galpha = TMath::ATan2(agxyz[1],agxyz[0]);
373	Double_t alpha = galpha;
374	galpha*=180/TMath::Pi();
375
376	if (galpha<45.) alpha = TMath::Pi()/8;
377	if (galpha>135.) alpha = TMath::Pi()*(7/8);
378	if (galpha>45 && galpha<135.) alpha = TMath::Pi()/2;
379
380	if (isOk) upper->Rotate(alpha);
381
382	if (isOk) isOk = upper->PropagateTo(kxAcorde,0);
383
384	TVectorD rgxyz(3);
385	AliExternalTrackParam param = (AliExternalTrackParam )upper;
386	param->GetXYZ(rgxyz.GetMatrixArray());
387	xAco = rgxyz[0];
388	yAco = rgxyz[1];
389	zAco = rgxyz[2];
390	fTracksToAcorde->Fill(xAco,zAco);
391	// Count how many tracks have been matched
392	nMatchedTracks++;
393	} else nSingleTracks++;
394	}
395
396	if ( (nMatchedTracks+nSingleTracks) <= nTracks ) muons=nMatchedTracks+nSingleTracks;
397	nTracksMatched->Fill(muons);
398
399	nMatchTracks = nSingleTracks + nMatchedTracks;
400
401	// Fills the tree
402
403	cosmicTree->Fill();
404
405	// Post output data.
406
407	PostData(0,cosmicTree);
408	PostData(1,histo);
409	}
410
411	//________________________________________________________________________
412	void AliAnalysisTaskAcorde::Terminate(Option_t *)
413	{
414
415	// AliDebug(1,"Do nothig in Terminate");
416
417	//nPart->Draw();
418	TCanvas *c1 = new TCanvas("c1","TPC-resolution for P reconstruction .... MRC");
419
420	c1->Divide(2,2);
421	c1->cd(1);
422	acordeHitsAll->Draw();
423	c1->cd(2)->SetLogy();
424	acordeMultiAll->Draw();
425	c1->cd(3)->SetLogy();
426	//nTracksMatched->Draw();
427	ntracks->Draw();
428	c1->cd(4);
429	fTracksToAcorde->Draw();
430
431
432
433
434	}