[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 fHistograms 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),
	fCosmicTree(0),
	fnTracks(0),
	fNMatchTracks(0),
	fkMinTPCclusters(30),
	fkMinTrackDist(1000.),
	fkMinCutDir(0.95),
  	fXAco(0),
  	fYAco(0),
  	fZAco(0),
	fTrigger(0),
	fActiveTriggerDetector(0),
	fNSingleTracks(0),
	fNMatchedTracks(0),
	fHisto(0),
 	fNTracks(0),
	fAcordeHitsAll(0),
	fAcordeMultiAll(0),
	fAcordeHitsTPC(0),
	fAcordeMultiTPC(0),
	fNTracksMatched(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 fHisto;
	delete fCosmicTree;
}
//________________________________________________________________________
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 fHistograms

	fHisto = new TList();

	fNTracksMatched = new TH1F("fNTracksMatched","fNTracksMatched from matching algorithm implemented",300,1,300);
	fNTracks = new TH1F("fNTracks","fnTracks distribution",500,1,500);
	fAcordeHitsAll = new TH1F("fAcordeHitsAll","Hits of ACORDE from ESD",60,-0.5,59.5);
	fAcordeHitsTPC = new TH1F("fAcordeHitsAllTPC","Hits of ACORDE from ESD (if fNTracks>0)",61,0,60);
	fAcordeMultiAll = new TH1F("fAcordeMultiAll","Multiplicity of ACORDE modules from ESD",60,-0.5,59.5);
	fAcordeMultiTPC = new TH1F("fAcordeMultiAllTPC","Multiplicity of ACORDE modules from ESD (id fNTracks>0)",61,0,60);
	fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);

	fHisto->Add(fNTracksMatched);
	fHisto->Add(fNTracks);
	fHisto->Add(fAcordeHitsAll);
	fHisto->Add(fAcordeHitsTPC);
	fHisto->Add(fAcordeMultiAll);
	fHisto->Add(fAcordeMultiTPC);
	fHisto->Add(fTracksToAcorde);

	// Create Tree branches
  	// Called just once

  	fCosmicTree = new TTree("fCosmicTree","fCosmicTreeMRC");
	fCosmicTree->Branch("fnTracks",&fnTracks,"fnTracks/I");
	fCosmicTree->Branch("fNMatchTracks",&fNMatchTracks,"fNMatchTracks/I");
	fCosmicTree->Branch("fXAco",&fXAco,"fXAco/F");
	fCosmicTree->Branch("fYAco",&fYAco,"fYAco/F");
	fCosmicTree->Branch("fZAco",&fZAco,"fZAco/F");
	fCosmicTree->Branch("fTrigger",&fTrigger,"fTrigger/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)) 
		{
			fAcordeHitsAll->Fill(imod);
			contMulti++;
		}

	}fAcordeMultiAll->Fill(contMulti);

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


	// Assingment of the type of Trigger to fCosmicTree

	fActiveTriggerDetector = fESD->GetFiredTriggerClasses();

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

	// Begins the matching analysis between tracks

	fnTracks = fESD->GetNumberOfTracks();
	if (fnTracks<=0) return;
	fNTracks->Fill(fnTracks);

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

		// Special case: fnTracks == 1
		if (fnTracks ==1)
		{
			AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
			if (!singleTrack || !singleTrack->GetOuterParam() || !singleTrack->GetInnerParam()) continue;
			if (singleTrack->GetTPCNcls()<fkMinTPCclusters) 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 
			fNMatchTracks=0;
			fNSingleTracks = fnTracks;
			muons = 1;

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

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

			Double_t trackdir0[3];
			track0->GetDirection(trackdir0);
			Float_t minDist = fkMinTrackDist;	

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

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

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

			for (Int_t jTrack=iTrack; jTrack<fnTracks; ++jTrack)
			{
				AliESDtrack *track1 = fESD->GetTrack(jTrack);
				if (!track1 || !track1->GetOuterParam() || !track1->GetInnerParam()) continue;
				Double_t trackdir1[3];
				track1->GetDirection(trackdir1);
				if (track1->GetTPCNcls()<fkMinTPCclusters) continue;
				Float_t direction = trackdir0[0]*trackdir1[0] + trackdir0[1]*trackdir1[1] + trackdir0[2]*trackdir1[2];
				if (TMath::Abs(direction)<fkMinCutDir) 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;


			}	
		} // fnTracks > 1
	

	} // Loop over all the tracks


	// Matching tracks
	
	for (Int_t itrack=0;itrack<fnTracks;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 kfXAcorde = 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(kfXAcorde,0);

			TVectorD rgxyz(3);
			AliExternalTrackParam *param = (AliExternalTrackParam *)upper;
			param->GetXYZ(rgxyz.GetMatrixArray());
			fXAco = rgxyz[0];
			fYAco = rgxyz[1];
			fZAco = rgxyz[2];
			fTracksToAcorde->Fill(fXAco,fZAco);
			// Count how many tracks have been matched
			fNMatchedTracks++;
		} else fNSingleTracks++;
	}

	if ( (fNMatchedTracks+fNSingleTracks) <= fnTracks ) muons=fNMatchedTracks+fNSingleTracks;
	fNTracksMatched->Fill(muons);

	fNMatchTracks = fNSingleTracks + fNMatchedTracks;
	
	// Fills the tree

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

	PostData(0,fCosmicTree);
	PostData(1,fHisto);
}      

//________________________________________________________________________
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);
	fAcordeHitsAll->Draw();
	c1->cd(2)->SetLogy();
	fAcordeMultiAll->Draw();
	c1->cd(3)->SetLogy();
	//fNTracksMatched->Draw();
	fNTracks->Draw();
	c1->cd(4);
	fTracksToAcorde->Draw();
	
	
}
Commit	Line	Data
	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 fHistograms 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	fCosmicTree(0),
	52	fnTracks(0),
	53	fNMatchTracks(0),
	54	fkMinTPCclusters(30),
	55	fkMinTrackDist(1000.),
	56	fkMinCutDir(0.95),
	57	fXAco(0),
	58	fYAco(0),
	59	fZAco(0),
	60	fTrigger(0),
	61	fActiveTriggerDetector(0),
	62	fNSingleTracks(0),
	63	fNMatchedTracks(0),
	64	fHisto(0),
	65	fNTracks(0),
	66	fAcordeHitsAll(0),
	67	fAcordeMultiAll(0),
	68	fAcordeHitsTPC(0),
	69	fAcordeMultiTPC(0),
	70	fNTracksMatched(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 fHisto;
	88	delete fCosmicTree;
	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 fHistograms
	116
	117	fHisto = new TList();
	118
	119	fNTracksMatched = new TH1F("fNTracksMatched","fNTracksMatched from matching algorithm implemented",300,1,300);
	120	fNTracks = new TH1F("fNTracks","fnTracks distribution",500,1,500);
	121	fAcordeHitsAll = new TH1F("fAcordeHitsAll","Hits of ACORDE from ESD",60,-0.5,59.5);
	122	fAcordeHitsTPC = new TH1F("fAcordeHitsAllTPC","Hits of ACORDE from ESD (if fNTracks>0)",61,0,60);
	123	fAcordeMultiAll = new TH1F("fAcordeMultiAll","Multiplicity of ACORDE modules from ESD",60,-0.5,59.5);
	124	fAcordeMultiTPC = new TH1F("fAcordeMultiAllTPC","Multiplicity of ACORDE modules from ESD (id fNTracks>0)",61,0,60);
	125	fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);
	126
	127	fHisto->Add(fNTracksMatched);
	128	fHisto->Add(fNTracks);
	129	fHisto->Add(fAcordeHitsAll);
	130	fHisto->Add(fAcordeHitsTPC);
	131	fHisto->Add(fAcordeMultiAll);
	132	fHisto->Add(fAcordeMultiTPC);
	133	fHisto->Add(fTracksToAcorde);
	134
	135	// Create Tree branches
	136	// Called just once
	137
	138	fCosmicTree = new TTree("fCosmicTree","fCosmicTreeMRC");
	139	fCosmicTree->Branch("fnTracks",&fnTracks,"fnTracks/I");
	140	fCosmicTree->Branch("fNMatchTracks",&fNMatchTracks,"fNMatchTracks/I");
	141	fCosmicTree->Branch("fXAco",&fXAco,"fXAco/F");
	142	fCosmicTree->Branch("fYAco",&fYAco,"fYAco/F");
	143	fCosmicTree->Branch("fZAco",&fZAco,"fZAco/F");
	144	fCosmicTree->Branch("fTrigger",&fTrigger,"fTrigger/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	fAcordeHitsAll->Fill(imod);
	178	contMulti++;
	179	}
	180
	181	}fAcordeMultiAll->Fill(contMulti);
	182
	183	for(Int_t imod=0;imod<60;imod++)
	184	{
	185	if (acordeESD->GetHitChannel(imod))
	186	{
	187	fAcordeHitsTPC->Fill(imod);
	188	contMultiTPC++;
	189	}
	190	}fAcordeMultiTPC->Fill(contMultiTPC);
	191
	192
	193	// Assingment of the type of Trigger to fCosmicTree
	194
	195	fActiveTriggerDetector = fESD->GetFiredTriggerClasses();
	196
	197	if (fActiveTriggerDetector.Contains("C0SCO")) fTrigger=0;
	198	if (fActiveTriggerDetector.Contains("C0SH2")) fTrigger=1;
	199	if (fActiveTriggerDetector.Contains("C0AMU")) fTrigger=2;
	200	if (fActiveTriggerDetector.Contains("C0LSR")) fTrigger=3;
	201	if (fActiveTriggerDetector.Contains("C0SCO1")) fTrigger=4;
	202	if (fActiveTriggerDetector.Contains("C0OBE")) fTrigger=5;
	203	if (fActiveTriggerDetector.Contains("C0PCO")) fTrigger=6;
	204	if (fActiveTriggerDetector.Contains("C0OB3")) fTrigger=7;
	205	if (fActiveTriggerDetector.Contains("C0OCP")) fTrigger=8;
	206	if (fActiveTriggerDetector.Contains("C0ASL")) fTrigger=9;
	207
	208
	209	// Begins the matching analysis between tracks
	210
	211	fnTracks = fESD->GetNumberOfTracks();
	212	if (fnTracks<=0) return;
	213	fNTracks->Fill(fnTracks);
	214
	215	fPair = new TArrayI(fnTracks);
	216	fNSingleTracks=0;
	217	fNMatchedTracks=0;
	218	Int_t muons = 0;
	219	TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
	220	for (Int_t iTrack=0; iTrack<fnTracks; iTrack++)
	221	{
	222
	223	// Special case: fnTracks == 1
	224	if (fnTracks ==1)
	225	{
	226	AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
	227	if (!singleTrack \|\| !singleTrack->GetOuterParam() \|\| !singleTrack->GetInnerParam()) continue;
	228	if (singleTrack->GetTPCNcls()<fkMinTPCclusters) 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	fNMatchTracks=0;
	235	fNSingleTracks = fnTracks;
	236	muons = 1;
	237
	238	} // fnTracks ==1
	239	else if (fnTracks>1) // fnTracks > 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 trackdir0[3];
	247	track0->GetDirection(trackdir0);
	248	Float_t minDist = fkMinTrackDist;
	249
	250	const AliExternalTrackParam * trackIn = track0->GetInnerParam();
	251	const AliExternalTrackParam * trackOut = track0->GetOuterParam();
	252
	253	if (!trackIn \|\| !trackOut) continue;
	254	if (fnTracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // Filter TPC-Laser
	255
	256	if (track0->GetTPCNcls()<fkMinTPCclusters) continue;
	257
	258	for (Int_t jTrack=iTrack; jTrack<fnTracks; ++jTrack)
	259	{
	260	AliESDtrack *track1 = fESD->GetTrack(jTrack);
	261	if (!track1 \|\| !track1->GetOuterParam() \|\| !track1->GetInnerParam()) continue;
	262	Double_t trackdir1[3];
	263	track1->GetDirection(trackdir1);
	264	if (track1->GetTPCNcls()<fkMinTPCclusters) continue;
	265	Float_t direction = trackdir0[0]trackdir1[0] + trackdir0[1]trackdir1[1] + trackdir0[2]*trackdir1[2];
	266	if (TMath::Abs(direction)<fkMinCutDir) 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	} // fnTracks > 1
	303
	304
	305	} // Loop over all the tracks
	306
	307
	308
	309
	310	// Matching tracks
	311
	312	for (Int_t itrack=0;itrack<fnTracks;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 kfXAcorde = 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(kfXAcorde,0);
	383
	384	TVectorD rgxyz(3);
	385	AliExternalTrackParam param = (AliExternalTrackParam )upper;
	386	param->GetXYZ(rgxyz.GetMatrixArray());
	387	fXAco = rgxyz[0];
	388	fYAco = rgxyz[1];
	389	fZAco = rgxyz[2];
	390	fTracksToAcorde->Fill(fXAco,fZAco);
	391	// Count how many tracks have been matched
	392	fNMatchedTracks++;
	393	} else fNSingleTracks++;
	394	}
	395
	396	if ( (fNMatchedTracks+fNSingleTracks) <= fnTracks ) muons=fNMatchedTracks+fNSingleTracks;
	397	fNTracksMatched->Fill(muons);
	398
	399	fNMatchTracks = fNSingleTracks + fNMatchedTracks;
	400
	401	// Fills the tree
	402
	403	fCosmicTree->Fill();
	404
	405	// Post output data.
	406
	407	PostData(0,fCosmicTree);
	408	PostData(1,fHisto);
	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	fAcordeHitsAll->Draw();
	423	c1->cd(2)->SetLogy();
	424	fAcordeMultiAll->Draw();
	425	c1->cd(3)->SetLogy();
	426	//fNTracksMatched->Draw();
	427	fNTracks->Draw();
	428	c1->cd(4);
	429	fTracksToAcorde->Draw();
	430
	431
	432
	433
	434	}