[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDtrackingEfficiencyCombined.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: AliTRDtrackingEfficiencyCombined.cxx 27496 2008-07-22 08:35:45Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Reconstruction QA                                                     //
//                                                                        //
//  Authors:                                                              //
//    Markus Fasel <M.Fasel@gsi.de>                                       //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <TObjArray.h>
#include <TProfile.h>
#include <TMath.h>
#include <TCanvas.h>
#include "TTreeStream.h"

#include "AliMagFMaps.h"
#include "AliTracker.h"
#include "AliTrackReference.h"
#include "AliAnalysisManager.h"

#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
#include "AliTRDtrackingEfficiencyCombined.h"

ClassImp(AliTRDtrackingEfficiencyCombined)

//_____________________________________________________________________________
AliTRDtrackingEfficiencyCombined::AliTRDtrackingEfficiencyCombined()
  :AliTRDrecoTask("TrackingEffMC", "Combined Tracking Efficiency")
{
  //
  // Default constructor
  //
}


//_____________________________________________________________________________
void AliTRDtrackingEfficiencyCombined::CreateOutputObjects(){
  //
  // Create output objects
  //

  OpenFile(0, "RECREATE");

  const Int_t nbins = 11;
  Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
  
  fContainer = new TObjArray();
  fContainer->Add(new TProfile("trEffComb", "Combined Tracking Efficiency", nbins, xbins));
  fContainer->Add(new TProfile("trContComb", "Combined Tracking Contamination", nbins, xbins));
}

//_____________________________________________________________________________
void AliTRDtrackingEfficiencyCombined::Exec(Option_t *){
  //
  // Do it
  //

	const Float_t kAlpha = 0.349065850;
	Int_t naccepted = 0, nrejected = 0, ndoublecounted = 0;
	Int_t labelsacc[10000];
	Int_t labelsrej[10000];
	Float_t momacc[10000];
	Float_t momrej[10000];
	TProfile *efficiency = (TProfile *)fContainer->At(0);
	TProfile *contamination = (TProfile *)fContainer->At(1);
	
	Int_t nTrackInfos = fTracks->GetEntriesFast();
	Double_t mom = 0;
	AliTRDtrackV1 *TRDtrack = 0x0;
	AliTRDtrackInfo *trkInf = 0x0;
	AliTrackReference *trackRef = 0x0;
	for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
		mom = 0.;
		trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
		if(!trkInf) continue;
		if((TRDtrack = trkInf->GetTRDtrack()) || trkInf->GetNumberOfClustersRefit()){
			// check if allready found by the tracker
			Bool_t found = kFALSE;
			for(Int_t il = 0; il < naccepted; il++){
				if(labelsacc[il] == trkInf->GetLabel()) found = kTRUE;
			}
			if(found){
				mom =  trackRef ? trackRef->P() : TRDtrack->P();
				contamination->Fill(mom, 1);
				ndoublecounted++;
				continue;
			}
			if(trkInf->GetNTrackRefs()){
				Int_t iref = 0;
				while(!(trackRef = trkInf->GetTrackRef(iref++)));
			}
			if(!trackRef) printf("Error: Track Reference missing for Track %d\n", trkInf->GetLabel());
			mom =  trackRef ? trackRef->P() : trkInf->GetOuterParam()->P();

//           Accept track
			if(fDebugLevel > 3)printf("Accept track\n");
			momacc[naccepted] = mom; 
			labelsacc[naccepted++] = trkInf->GetLabel();
/*			printf("Reconstructed: event %3d Tracks: MC[%d] ESD[%d] NRefs[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), trkInf->GetLabel(),  trkInf->GetTrackId(), trkInf->GetNTrackRefs());*/
    } else{
			if(fDebugLevel>10) printf("Analysing Track References\n");
      // Check if track is findable
			Float_t xmin = 10000.0, xmax = 0.0; 
			Float_t ymin = 0.0, ymax = 0.0;
			Float_t zmin = 0.0, zmax = 0.0;
			Float_t lastx = 0.0, x = 0.0;
			Int_t nLayers = 0;
/*			trackRef = trkInf->GetTrackRef(0);*/
/*			xmin = trackRef->LocalX(); xmax = trackRef->LocalX();
      ymin = trackRef->LocalY(); ymax = trackRef->LocalY();
      mom = trackRef->P();*/
      Int_t *sector = new Int_t[trkInf->GetNTrackRefs()];
      for(Int_t itr = 0; itr < trkInf->GetNTrackRefs(); itr++){
        trackRef = trkInf->GetTrackRef(itr);
				if(fDebugLevel>10) printf("%d. x[%f], y[%f], z[%f]\n", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z());
				x = trackRef->LocalX(); 
				if(x < 250. || x > 370.) continue;	// Be Sure that we are inside TRD
        sector[itr] = Int_t(trackRef->Alpha()/kAlpha);
        if(x < xmin){
          xmin = trackRef->LocalX();
          ymin = trackRef->LocalY();
          zmin = trackRef->Z();
          mom = trackRef->P();
        }
        if(x > xmax){
          xmax = trackRef->LocalX();
          ymax = trackRef->LocalY();
          zmax = trackRef->Z();
        }
				if(itr > 0){
					Float_t dist = TMath::Abs(x - lastx);
					if(fDebugLevel>10) printf("x = %f, lastx = %f, dist = %f\n", x, lastx, dist);
					if(TMath::Abs(dist - 3.7) < 0.1) nLayers++; 	// ref(i+1) has to be larger than ref(i)
				}
				lastx = x;
      }
      // Apply cuts
      Bool_t findable = kTRUE;
      if(trkInf->GetNTrackRefs() > 2 && xmax > xmin){
        if(mom < 0.55) findable = kFALSE;									// momentum cut at 0.6
 				Double_t yangle = (ymax -ymin)/(xmax - xmin);
				Double_t zangle = (zmax -zmin)/(xmax - xmin);
				if(fDebugLevel>10) printf("track: y-Angle = %f, z-Angle = %f\n", yangle, zangle);
				if(fDebugLevel>10) printf("nLayers = %d\n", nLayers);
				if(TMath::ATan(TMath::Abs(yangle)) > 45.) findable = kFALSE;
				if(TMath::ATan(TMath::Abs(zangle)) > 45.) findable = kFALSE;
				if(nLayers < 4) findable = kFALSE;
        if(!trkInf->IsPrimary()) findable = kFALSE;
        Bool_t samesec = kTRUE;
        for(Int_t iref = 1; iref < trkInf->GetNTrackRefs(); iref++)
          if(sector[iref] != sector[0]) samesec = kFALSE;
        if(!samesec) findable = kFALSE;		// Discard all tracks which are passing more than one sector
        if(fDebugLevel){
          Double_t trackAngle = TMath::ATan(yangle);
          Bool_t primary = trkInf->IsPrimary();
          (*fDebugStream) << "NotFoundTrack"
            << "Momentum=" 	<< mom
            << "trackAngle="<< trackAngle
            << "NLayers="	<< nLayers
            << "Primary="	<< primary
            << "\n";
        }
      }
      else
        findable = kFALSE;
      delete[] sector;
      if(findable){
				momrej[nrejected] = mom;
				labelsrej[nrejected++] = trkInf->GetLabel();
/*			  printf("Not Reconstructed: event %3d Tracks: MC[%d] ESD[%d] NRefs[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), trkInf->GetLabel(),  trkInf->GetTrackId(), trkInf->GetNTrackRefs());*/
      }
    }
  }
	for(Int_t itk = 0; itk < naccepted; itk++){
		if(fDebugLevel > 2)printf("Accepted MC track: %d\n", labelsacc[itk]);
		efficiency->Fill(momacc[itk], 1);
		contamination->Fill(momacc[itk], 0);
	}
	Int_t nall = naccepted;
	for(Int_t imis = 0; imis < nrejected; imis++){
		Bool_t found = kFALSE;
		for(Int_t ifound = 0; ifound < naccepted; ifound++){
			if(labelsacc[ifound] == labelsrej[imis]){
				found = kTRUE;
				break;
			}
		}
		if(!found){
			efficiency->Fill(momrej[imis], 0);
			contamination->Fill(momrej[imis], 0);
			if(fDebugLevel > 2)printf("Rejected MC track: %d\n", labelsrej[imis]);
			nall++;
		}
	}
  //if(fDebugLevel>=1)
  printf("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall, naccepted, 1.E2*Float_t(naccepted)/Float_t(nall));
  printf("%3d Tracks: ALL[%3d] DoubleCounted[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall + ndoublecounted, ndoublecounted, 1.E2*Float_t(ndoublecounted)/Float_t(nall + ndoublecounted));

  PostData(0, fContainer);
}

//_____________________________________________________________________________
void AliTRDtrackingEfficiencyCombined::Terminate(Option_t *)
{
  //
  // Termination
  //

  if(fDebugStream){ 
    delete fDebugStream;
    fDebugStream = 0x0;
    fDebugLevel = 0;
  }

  fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
  if (!fContainer) {
    Printf("ERROR: list not available");
    return;
  }

  /*TProfile *hEff = (TProfile*)fContainer->At(0);
  TProfile *hEffCont = (TProfile*)fContainer->At(1);
  Printf("Eff[%p] EffCont[%p]\n", (void*)hEff, (void*)hEffCont);


  TCanvas *c2 = new TCanvas("c2","",800,400);
  c2->Divide(2,1);

  c2->cd(1);
  hEff->DrawCopy("e1");
  c2->cd(2);
  hEffCont->DrawCopy("e1");*/
}
Commit	Line	Data
814ecea4	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id: AliTRDtrackingEfficiencyCombined.cxx 27496 2008-07-22 08:35:45Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Reconstruction QA //
	21	// //
	22	// Authors: //
	23	// Markus Fasel <M.Fasel@gsi.de> //
	24	// //
	25	////////////////////////////////////////////////////////////////////////////
	26
	27	#include <TObjArray.h>
	28	#include <TProfile.h>
	29	#include <TMath.h>
	30	#include <TCanvas.h>
3d86166d	31	#include "TTreeStream.h"
814ecea4	32
	33	#include "AliMagFMaps.h"
	34	#include "AliTracker.h"
814ecea4	35	#include "AliTrackReference.h"
814ecea4	36	#include "AliAnalysisManager.h"
3d86166d	37
	38	#include "AliTRDseedV1.h"
	39	#include "AliTRDtrackV1.h"
814ecea4	40	#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
	41	#include "AliTRDtrackingEfficiencyCombined.h"
	42
	43	ClassImp(AliTRDtrackingEfficiencyCombined)
	44
	45	//_____________________________________________________________________________
3d86166d	46	AliTRDtrackingEfficiencyCombined::AliTRDtrackingEfficiencyCombined()
3d86166d	47	:AliTRDrecoTask("TrackingEffMC", "Combined Tracking Efficiency")
814ecea4	48	{
	49	//
	50	// Default constructor
	51	//
ed383798	52	}
ed383798	53
814ecea4	54
	55	//_____________________________________________________________________________
	56	void AliTRDtrackingEfficiencyCombined::CreateOutputObjects(){
	57	//
	58	// Create output objects
	59	//
	60
	61	OpenFile(0, "RECREATE");
	62
	63	const Int_t nbins = 11;
	64	Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
	65
3d86166d	66	fContainer = new TObjArray();
	67	fContainer->Add(new TProfile("trEffComb", "Combined Tracking Efficiency", nbins, xbins));
	68	fContainer->Add(new TProfile("trContComb", "Combined Tracking Contamination", nbins, xbins));
814ecea4	69	}
	70
	71	//_____________________________________________________________________________
	72	void AliTRDtrackingEfficiencyCombined::Exec(Option_t *){
	73	//
	74	// Do it
	75	//
	76
	77	const Float_t kAlpha = 0.349065850;
	78	Int_t naccepted = 0, nrejected = 0, ndoublecounted = 0;
	79	Int_t labelsacc[10000];
	80	Int_t labelsrej[10000];
	81	Float_t momacc[10000];
	82	Float_t momrej[10000];
3d86166d	83	TProfile efficiency = (TProfile )fContainer->At(0);
3d86166d	84	TProfile contamination = (TProfile )fContainer->At(1);
814ecea4	85
3d86166d	86	Int_t nTrackInfos = fTracks->GetEntriesFast();
814ecea4	87	Double_t mom = 0;
	88	AliTRDtrackV1 *TRDtrack = 0x0;
	89	AliTRDtrackInfo *trkInf = 0x0;
	90	AliTrackReference *trackRef = 0x0;
	91	for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
	92	mom = 0.;
3d86166d	93	trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
814ecea4	94	if(!trkInf) continue;
	95	if((TRDtrack = trkInf->GetTRDtrack()) \|\| trkInf->GetNumberOfClustersRefit()){
	96	// check if allready found by the tracker
	97	Bool_t found = kFALSE;
	98	for(Int_t il = 0; il < naccepted; il++){
	99	if(labelsacc[il] == trkInf->GetLabel()) found = kTRUE;
	100	}
	101	if(found){
	102	mom = trackRef ? trackRef->P() : TRDtrack->P();
	103	contamination->Fill(mom, 1);
	104	ndoublecounted++;
	105	continue;
	106	}
	107	if(trkInf->GetNTrackRefs()){
	108	Int_t iref = 0;
	109	while(!(trackRef = trkInf->GetTrackRef(iref++)));
	110	}
158754f9	111	if(!trackRef) printf("Error: Track Reference missing for Track %d\n", trkInf->GetLabel());
b0a48c4d	112	mom = trackRef ? trackRef->P() : trkInf->GetOuterParam()->P();
b0a48c4d	113
814ecea4	114	// Accept track
	115	if(fDebugLevel > 3)printf("Accept track\n");
	116	momacc[naccepted] = mom;
	117	labelsacc[naccepted++] = trkInf->GetLabel();
	118	/* printf("Reconstructed: event %3d Tracks: MC[%d] ESD[%d] NRefs[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), trkInf->GetLabel(), trkInf->GetTrackId(), trkInf->GetNTrackRefs());*/
	119	} else{
	120	if(fDebugLevel>10) printf("Analysing Track References\n");
	121	// Check if track is findable
	122	Float_t xmin = 10000.0, xmax = 0.0;
	123	Float_t ymin = 0.0, ymax = 0.0;
	124	Float_t zmin = 0.0, zmax = 0.0;
	125	Float_t lastx = 0.0, x = 0.0;
	126	Int_t nLayers = 0;
	127	/* trackRef = trkInf->GetTrackRef(0);*/
	128	/* xmin = trackRef->LocalX(); xmax = trackRef->LocalX();
	129	ymin = trackRef->LocalY(); ymax = trackRef->LocalY();
	130	mom = trackRef->P();*/
	131	Int_t *sector = new Int_t[trkInf->GetNTrackRefs()];
	132	for(Int_t itr = 0; itr < trkInf->GetNTrackRefs(); itr++){
	133	trackRef = trkInf->GetTrackRef(itr);
	134	if(fDebugLevel>10) printf("%d. x[%f], y[%f], z[%f]\n", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z());
	135	x = trackRef->LocalX();
	136	if(x < 250. \|\| x > 370.) continue; // Be Sure that we are inside TRD
	137	sector[itr] = Int_t(trackRef->Alpha()/kAlpha);
	138	if(x < xmin){
	139	xmin = trackRef->LocalX();
	140	ymin = trackRef->LocalY();
	141	zmin = trackRef->Z();
	142	mom = trackRef->P();
	143	}
	144	if(x > xmax){
	145	xmax = trackRef->LocalX();
	146	ymax = trackRef->LocalY();
	147	zmax = trackRef->Z();
	148	}
	149	if(itr > 0){
	150	Float_t dist = TMath::Abs(x - lastx);
	151	if(fDebugLevel>10) printf("x = %f, lastx = %f, dist = %f\n", x, lastx, dist);
	152	if(TMath::Abs(dist - 3.7) < 0.1) nLayers++; // ref(i+1) has to be larger than ref(i)
	153	}
	154	lastx = x;
	155	}
	156	// Apply cuts
	157	Bool_t findable = kTRUE;
	158	if(trkInf->GetNTrackRefs() > 2 && xmax > xmin){
	159	if(mom < 0.55) findable = kFALSE; // momentum cut at 0.6
	160	Double_t yangle = (ymax -ymin)/(xmax - xmin);
	161	Double_t zangle = (zmax -zmin)/(xmax - xmin);
	162	if(fDebugLevel>10) printf("track: y-Angle = %f, z-Angle = %f\n", yangle, zangle);
	163	if(fDebugLevel>10) printf("nLayers = %d\n", nLayers);
	164	if(TMath::ATan(TMath::Abs(yangle)) > 45.) findable = kFALSE;
	165	if(TMath::ATan(TMath::Abs(zangle)) > 45.) findable = kFALSE;
	166	if(nLayers < 4) findable = kFALSE;
	167	if(!trkInf->IsPrimary()) findable = kFALSE;
	168	Bool_t samesec = kTRUE;
	169	for(Int_t iref = 1; iref < trkInf->GetNTrackRefs(); iref++)
	170	if(sector[iref] != sector[0]) samesec = kFALSE;
	171	if(!samesec) findable = kFALSE; // Discard all tracks which are passing more than one sector
	172	if(fDebugLevel){
	173	Double_t trackAngle = TMath::ATan(yangle);
	174	Bool_t primary = trkInf->IsPrimary();
	175	(*fDebugStream) << "NotFoundTrack"
	176	<< "Momentum=" << mom
	177	<< "trackAngle="<< trackAngle
178	<< "NLayers=" << nLayers
179	<< "Primary=" << primary
180	<< "\n";
181	}
182	}
183	else
184	findable = kFALSE;
185	delete[] sector;
186	if(findable){
187	momrej[nrejected] = mom;
188	labelsrej[nrejected++] = trkInf->GetLabel();
189	/* printf("Not Reconstructed: event %3d Tracks: MC[%d] ESD[%d] NRefs[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), trkInf->GetLabel(), trkInf->GetTrackId(), trkInf->GetNTrackRefs());*/
190	}
191	}
192	}
193	for(Int_t itk = 0; itk < naccepted; itk++){
194	if(fDebugLevel > 2)printf("Accepted MC track: %d\n", labelsacc[itk]);
195	efficiency->Fill(momacc[itk], 1);
196	contamination->Fill(momacc[itk], 0);
197	}
198	Int_t nall = naccepted;
199	for(Int_t imis = 0; imis < nrejected; imis++){
200	Bool_t found = kFALSE;
201	for(Int_t ifound = 0; ifound < naccepted; ifound++){
202	if(labelsacc[ifound] == labelsrej[imis]){
203	found = kTRUE;
204	break;
205	}
206	}
207	if(!found){
208	efficiency->Fill(momrej[imis], 0);
209	contamination->Fill(momrej[imis], 0);
210	if(fDebugLevel > 2)printf("Rejected MC track: %d\n", labelsrej[imis]);
211	nall++;
212	}
213	}
214	//if(fDebugLevel>=1)
215	printf("%3d Tracks: MC[%3d] TRD[%3d \| %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall, naccepted, 1.E2*Float_t(naccepted)/Float_t(nall));
216	printf("%3d Tracks: ALL[%3d] DoubleCounted[%3d \| %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall + ndoublecounted, ndoublecounted, 1.E2*Float_t(ndoublecounted)/Float_t(nall + ndoublecounted));
217
3d86166d	218	PostData(0, fContainer);
814ecea4	219	}
	220
	221	//_____________________________________________________________________________
	222	void AliTRDtrackingEfficiencyCombined::Terminate(Option_t *)
	223	{
	224	//
	225	// Termination
	226	//
	227
3d86166d	228	if(fDebugStream){
	229	delete fDebugStream;
	230	fDebugStream = 0x0;
	231	fDebugLevel = 0;
	232	}
	233
	234	fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
	235	if (!fContainer) {
814ecea4	236	Printf("ERROR: list not available");
	237	return;
	238	}
	239
3d86166d	240	/TProfile hEff = (TProfile*)fContainer->At(0);
3d86166d	241	TProfile hEffCont = (TProfile)fContainer->At(1);
814ecea4	242	Printf("Eff[%p] EffCont[%p]\n", (void)hEff, (void)hEffCont);
	243
	244
	245	TCanvas *c2 = new TCanvas("c2","",800,400);
	246	c2->Divide(2,1);
	247
	248	c2->cd(1);
	249	hEff->DrawCopy("e1");
	250	c2->cd(2);
	251	hEffCont->DrawCopy("e1");*/
814ecea4	252	}