[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
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	//
6399c648	11	// Create some fHistograms and one tree with the information
b5f2577c	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),
6399c648	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),
b5f2577c	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;
6399c648	87	delete fHisto;
6399c648	88	delete fCosmicTree;
b5f2577c	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
6399c648	115	// creates one TList with some fHistograms
b5f2577c	116
6399c648	117	fHisto = new TList();
b5f2577c	118
6399c648	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);
b5f2577c	125	fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);
b5f2577c	126
6399c648	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);
b5f2577c	134
	135	// Create Tree branches
	136	// Called just once
	137
6399c648	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");
b5f2577c	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	{
6399c648	177	fAcordeHitsAll->Fill(imod);
b5f2577c	178	contMulti++;
	179	}
	180
6399c648	181	}fAcordeMultiAll->Fill(contMulti);
b5f2577c	182
	183	for(Int_t imod=0;imod<60;imod++)
	184	{
	185	if (acordeESD->GetHitChannel(imod))
	186	{
6399c648	187	fAcordeHitsTPC->Fill(imod);
b5f2577c	188	contMultiTPC++;
b5f2577c	189	}
6399c648	190	}fAcordeMultiTPC->Fill(contMultiTPC);
b5f2577c	191
b5f2577c	192
6399c648	193	// Assingment of the type of Trigger to fCosmicTree
b5f2577c	194
6399c648	195	fActiveTriggerDetector = fESD->GetFiredTriggerClasses();
b5f2577c	196
6399c648	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;
b5f2577c	207
	208
	209	// Begins the matching analysis between tracks
	210
6399c648	211	fnTracks = fESD->GetNumberOfTracks();
	212	if (fnTracks<=0) return;
	213	fNTracks->Fill(fnTracks);
b5f2577c	214
6399c648	215	fPair = new TArrayI(fnTracks);
	216	fNSingleTracks=0;
	217	fNMatchedTracks=0;
b5f2577c	218	Int_t muons = 0;
b5f2577c	219	TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
6399c648	220	for (Int_t iTrack=0; iTrack<fnTracks; iTrack++)
b5f2577c	221	{
b5f2577c	222
6399c648	223	// Special case: fnTracks == 1
6399c648	224	if (fnTracks ==1)
b5f2577c	225	{
	226	AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
	227	if (!singleTrack \|\| !singleTrack->GetOuterParam() \|\| !singleTrack->GetInnerParam()) continue;
6399c648	228	if (singleTrack->GetTPCNcls()<fkMinTPCclusters) continue;
b5f2577c	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
6399c648	234	fNMatchTracks=0;
6399c648	235	fNSingleTracks = fnTracks;
b5f2577c	236	muons = 1;
b5f2577c	237
6399c648	238	} // fnTracks ==1
6399c648	239	else if (fnTracks>1) // fnTracks > 1
b5f2577c	240	{
	241
	242	AliESDtrack *track0 = fESD->GetTrack(iTrack);
	243	(*fPair)[iTrack]=-1;
	244	if (!track0 \|\| !track0->GetOuterParam() \|\| !track0->GetInnerParam()) continue;
	245
6399c648	246	Double_t trackdir0[3];
	247	track0->GetDirection(trackdir0);
	248	Float_t minDist = fkMinTrackDist;
b5f2577c	249
	250	const AliExternalTrackParam * trackIn = track0->GetInnerParam();
	251	const AliExternalTrackParam * trackOut = track0->GetOuterParam();
	252
	253	if (!trackIn \|\| !trackOut) continue;
6399c648	254	if (fnTracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // Filter TPC-Laser
b5f2577c	255
6399c648	256	if (track0->GetTPCNcls()<fkMinTPCclusters) continue;
b5f2577c	257
6399c648	258	for (Int_t jTrack=iTrack; jTrack<fnTracks; ++jTrack)
b5f2577c	259	{
	260	AliESDtrack *track1 = fESD->GetTrack(jTrack);
	261	if (!track1 \|\| !track1->GetOuterParam() \|\| !track1->GetInnerParam()) continue;
6399c648	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
b5f2577c	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	}
6399c648	302	} // fnTracks > 1
b5f2577c	303
	304
	305	} // Loop over all the tracks
	306
	307
	308
	309
	310	// Matching tracks
	311
6399c648	312	for (Int_t itrack=0;itrack<fnTracks;itrack++)
b5f2577c	313	{
6399c648	314	AliESDtrack *upTrack = 0;
6399c648	315	AliESDtrack *downTrack = 0;
b5f2577c	316
6399c648	317	AliESDtrack *track0 = fESD->GetTrack(itrack);
6399c648	318	if (!track0 \|\| !track0->GetOuterParam() \|\| !track0->GetInnerParam() ) continue;
b5f2577c	319
b5f2577c	320	Double_t mxyz[3];
6399c648	321	track0->GetOuterParam()->GetXYZ(mxyz);
b5f2577c	322
6399c648	323	upTrack = track0;
b5f2577c	324
	325
	326	if ((*fPair)[itrack]>0)
	327	{
6399c648	328	AliESDtrack track1 = fESD->GetTrack((fPair)[itrack]);
6399c648	329	if (!track1 \|\| !track1->GetOuterParam() \|\| !track1->GetInnerParam()) continue;
b5f2577c	330	Float_t dvertex0[3];
	331	Float_t dvertex1[3];
	332
6399c648	333	dvertex0[0] = track0->Xv();
	334	dvertex0[1] = track0->Yv();
	335	dvertex0[2] = track0->Zv();
b5f2577c	336
6399c648	337	dvertex1[0] = track1->Xv();
	338	dvertex1[1] = track1->Yv();
	339	dvertex1[2] = track1->Zv();
b5f2577c	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
6399c648	344	if (track1->GetOuterParam())
b5f2577c	345	{
b5f2577c	346	Double_t nxyz[3];
6399c648	347	track1->GetOuterParam()->GetXYZ(nxyz);
b5f2577c	348	if (nxyz[1]>mxyz[1])
b5f2577c	349	{
6399c648	350	upTrack = track1;
6399c648	351	}else downTrack = track1;
b5f2577c	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
6399c648	360	const Double_t kfXAcorde = 850.; // radios of "ACORDE" above the magnet
b5f2577c	361
	362	Double_t agxyz[3];
	363
	364
6399c648	365	AliExternalTrackParam upper = (AliExternalTrackParam )(upTrack->GetOuterParam()->Clone());
b5f2577c	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
6399c648	382	if (isOk) isOk = upper->PropagateTo(kfXAcorde,0);
b5f2577c	383
	384	TVectorD rgxyz(3);
	385	AliExternalTrackParam param = (AliExternalTrackParam )upper;
	386	param->GetXYZ(rgxyz.GetMatrixArray());
6399c648	387	fXAco = rgxyz[0];
	388	fYAco = rgxyz[1];
	389	fZAco = rgxyz[2];
	390	fTracksToAcorde->Fill(fXAco,fZAco);
b5f2577c	391	// Count how many tracks have been matched
6399c648	392	fNMatchedTracks++;
6399c648	393	} else fNSingleTracks++;
b5f2577c	394	}
b5f2577c	395
6399c648	396	if ( (fNMatchedTracks+fNSingleTracks) <= fnTracks ) muons=fNMatchedTracks+fNSingleTracks;
6399c648	397	fNTracksMatched->Fill(muons);
b5f2577c	398
6399c648	399	fNMatchTracks = fNSingleTracks + fNMatchedTracks;
b5f2577c	400
	401	// Fills the tree
	402
6399c648	403	fCosmicTree->Fill();
b5f2577c	404
	405	// Post output data.
	406
6399c648	407	PostData(0,fCosmicTree);
6399c648	408	PostData(1,fHisto);
b5f2577c	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);
6399c648	422	fAcordeHitsAll->Draw();
b5f2577c	423	c1->cd(2)->SetLogy();
6399c648	424	fAcordeMultiAll->Draw();
b5f2577c	425	c1->cd(3)->SetLogy();
6399c648	426	//fNTracksMatched->Draw();
6399c648	427	fNTracks->Draw();
b5f2577c	428	c1->cd(4);
	429	fTracksToAcorde->Draw();
	430
	431
	432
	433
	434	}