[u/mrichter/AliRoot.git] / PWGPP / comparison / AliITSComparisonTask.cxx

//-------------------------------------------------------------------------
//
// This is the PROOF-enabled version of ITS/AliITSComparisonV2.C macro.
// Origin:  Andrei.Zalite@cern.ch
//
//-------------------------------------------------------------------------

#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TList.h"
#include "TClonesArray.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLine.h"
#include "TText.h"
#include "TFile.h"

#include "AliLog.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDtrack.h"
#include "AliTrackReference.h"
#include "AliMCComparisonTrack.h"

#include "AliITSComparisonTask.h"

ClassImp(AliITSComparisonTask)

extern TStyle *gStyle;

AliITSComparisonTask::AliITSComparisonTask()
  : AliAnalysisTaskSE("AliITSComaprisonTask"),
    fListOfHistos(0),
    fGood(0),
    fFound(0),
    fFake(0),
    fP(0),
    fL(0),
    fPt(0),
    fTip(0),
    fE(0),
    fEp(0),
    fGoodPhi(0),
    fFoundPhi(0),
    fLip(0)
{
  // Default constructor
  AliInfo("Default constructor AliITSComparisonTask");
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #1 TList
  DefineOutput(1, TList::Class());
}


AliITSComparisonTask::AliITSComparisonTask(const char* name)
  : AliAnalysisTaskSE(name),
    fListOfHistos(0),
    fGood(0),
    fFound(0),
    fFake(0),
    fP(0),
    fL(0),
    fPt(0),
    fTip(0),
    fE(0),
    fEp(0),
    fGoodPhi(0),
    fFoundPhi(0),
    fLip(0)
{
  // Constructor
  AliInfo("Constructor AliITSComparisonTask");
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #1 TList
  DefineOutput(1, TList::Class());
}


void AliITSComparisonTask::UserCreateOutputObjects() 
{
  // Create histograms
  // Called once
  AliInfo("AliITSComparisonTask::UserCreateOutputObjects");
  // Create output container
  fListOfHistos = new TList();
  
  fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
  fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
  fFake = new TH1F("fFake",  "Pt for fake tracks", 34, 0.2, 7.0);
  fP = new TH1F("fP", "PHI resolution", 50, -20., 20.);
  fL = new TH1F("fL", "LAMBDA resolution", 50, -20., 20.);
  fPt = new TH1F("fPt", "Relative Pt resolution", 30, -10., 10.);
  fTip = new TH1F("fTip", "Transverse impact parameter", 40, -800., 800.);
  fE = new TH1F("fE", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
  fEp = new TH2F("fEp", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
  fGoodPhi = new TH1F("fGoodPhi", "Phi for good tracks", 90, 0., 2.*TMath::Pi());
  fFoundPhi = new TH1F("fFoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
  fLip = new TH1F("fLip","Longitudinal impact parameter", 40, -800., 800.);
  
  
  fListOfHistos->Add(fGood);
  fListOfHistos->Add(fFound);
  fListOfHistos->Add(fFake);
  fListOfHistos->Add(fP);
  fListOfHistos->Add(fL);
  fListOfHistos->Add(fPt);
  fListOfHistos->Add(fTip);
  fListOfHistos->Add(fE);
  fListOfHistos->Add(fEp);
  fListOfHistos->Add(fGoodPhi);
  fListOfHistos->Add(fFoundPhi);
  fListOfHistos->Add(fLip);
}


void AliITSComparisonTask::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  
    // MC information
  AliMCEvent* mcEvent = MCEvent();
  if (!mcEvent) {
    Printf("ERROR: Could not retrieve MC event");
    return;
  }
   
  Int_t nt = 0;
    
  TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
    
  // Loop over all MC tracks and select "good" tracks
  for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
    AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
    if (!track) {
      Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
      continue;
    }
    
    // Track selection
    if (track->Pt() < 0.2) continue;
    if (track->Pt() > 7.0) continue;    
    if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
    
    Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
    if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
    if (TMath::Abs(vz) > 50.) continue; 
    
    // Loop over Track References
    Bool_t labelTPC = kFALSE;
    AliTrackReference* trackRef = 0;
    UInt_t iTSLayerMap = 0;
    for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
      trackRef = track->GetTrackReference(iTrackRef);
      if(trackRef) {
	Int_t detectorId = trackRef->DetectorId();
	if (detectorId == AliTrackReference::kITS) {
	  Float_t radius = trackRef->R();
	  if (radius > 2.5 && radius < 5.5)
	    iTSLayerMap |= 1;
	  else if (radius > 5.5 && radius < 8.5)
	    iTSLayerMap |= 1 << 1;
	  else if (radius > 13.5 && radius < 16.5)
	    iTSLayerMap |= 1 << 2;
	  else if (radius > 22. && radius < 26.)
	    iTSLayerMap |= 1 << 3;
	  else if (radius > 36. && radius < 41.)
	    iTSLayerMap |= 1 << 4;
	  else if (radius > 42. && radius < 46.)
	    iTSLayerMap |= 1 << 5;
	  else {
	    Printf("Wrong radius %f ", radius);
	    return;
	  }

	}
	if (detectorId == AliTrackReference::kTPC) {	    
	  labelTPC = kTRUE;
	  break;
	}
      }      
    }    // track references loop   
    
    // Skip tracks that passed not all ITS layers 
    if (iTSLayerMap != 0x3F) continue;
    
    // "Good" tracks
    if (labelTPC) {
      AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
      ref->SetLabel(iTracks);
      TParticle* particle = track->Particle();
      Int_t pdg = particle->GetPdgCode();
      ref->SetPDG(pdg); 
      ref->SetPz(track->Pz());
      Float_t pt = track->Pt();
      ref->SetPt(pt);
      fGood->Fill(pt);
      Float_t phig = track->Phi();
      ref->SetPhi(phig);
      fGoodPhi->Fill(phig);
      nt++;  
    }  
  }    // track loop 
  
  
  // ESD information  
  AliVEvent* event = InputEvent();
  if (!event) {
    Printf("ERROR: Could not retrieve event");
    return;
  }
    
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
  
  Bool_t iFound;
  Int_t nfound = 0;
  Int_t nfake = 0;
  Int_t nlost = 0;
  
  Int_t mcGoods = refs->GetEntriesFast();
  
  // Loop over all "good" MC tracks
  for (Int_t k = 0; k < mcGoods; k++) {
    AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k); 
    Int_t mcLabel = ref->GetLabel();
    Float_t ptg = ref->GetPt(); 
    Float_t phiG = ref->GetPhi();
    iFound = kFALSE;
    for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
      AliESDtrack* track = esd->GetTrack(iTrack);  
      if (!track) {
        Printf("ERROR: Could not receive track %d", iTrack);
        continue;
      }
      
      //      if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;
      if (! track->IsOn(AliESDtrack::kITSrefit)) continue;
      
      Int_t label =  track->GetLabel();
      
      if (mcLabel == TMath::Abs(label)) {	  
        if (mcLabel == label) {
	  nfound++;
	  fFound->Fill(ptg);
	  fFoundPhi->Fill(phiG);
	} 
	else {
	  nfake++;
	  fFake->Fill(ptg);
	}
	iFound = kTRUE;

	Float_t phi = track->Phi();
	Double_t pt = track->Pt();
	
        Float_t lam = TMath::ATan2(track->Pz(),pt); 
        Float_t pt1 = 1/pt;
	
	Float_t phig = ref->GetPhi();
        fP->Fill((phi - phig)*1000.);   
	
	Float_t lamg = TMath::ATan2(ref->GetPz(),ptg);
        fL->Fill((lam - lamg)*1000.);
	
        fPt->Fill((pt1 - 1/ptg)/(1/ptg)*100.);  
	
        Float_t d,z; track->GetImpactParameters(d,z);
        fTip->Fill(10000.*d);
        fLip->Fill(10000.*z);
	
	Float_t mom = pt/TMath::Cos(lam);
        Float_t dedx = track->GetTPCsignal();
        fEp->Fill(mom, dedx, 1.);
	
	Int_t pdg = ref->GetPDG();
        if (TMath::Abs(pdg )== 211) { // pions
	  if (mom > 0.4 && mom < 0.5) {
            fE->Fill(dedx, 1.);
          }
	}   
	   
        break; 
      }
    }  
    if (!iFound) {
      nlost++;
    } 
  }
     
  Printf(" Results: " );
  Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
  
  refs->Clear();
  
  // Post output data.
  PostData(1, fListOfHistos);
}      

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

  // Draw result to the screen
  // Called once at the end of the query 
  fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
  if (!fListOfHistos) {
    Printf("ERROR: fListOfHistos not available");
    return;
  }

  fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
  fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
  fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));  
  fP = dynamic_cast<TH1F*>(fListOfHistos->At(3));
  fL = dynamic_cast<TH1F*>(fListOfHistos->At(4));
  fPt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
  fTip = dynamic_cast<TH1F*>(fListOfHistos->At(6));
  fE = dynamic_cast<TH1F*>(fListOfHistos->At(7));
  fEp = dynamic_cast<TH2F*>(fListOfHistos->At(8));
  fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
  fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
  fLip = dynamic_cast<TH1F*>(fListOfHistos->At(11));
  
  gStyle->SetOptStat(111110);
  gStyle->SetOptFit(1);
  
  TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
  
  Int_t minc = 33;
  
  TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
  p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
  fP->SetFillColor(4); fP->SetXTitle("(mrad)");
  if (fP->GetEntries() < minc) fP->Draw(); else fP->Fit("gaus"); c1->cd();
  
  TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1., 0.6); p2->Draw();
  p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
  fL->SetFillColor(4); fL->SetXTitle("(mrad)");
  if (fL->GetEntries() < minc) fL->Draw(); else fL->Fit("gaus"); c1->cd();
  
  TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
  p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
  fPt->SetFillColor(2); fPt->SetXTitle("(%)");
  if (fPt->GetEntries() < minc) fPt->Draw(); else fPt->Fit("gaus"); c1->cd();
  
  TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
  p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
  fTip->SetFillColor(6); fTip->SetXTitle("(micron)");
  if (fTip->GetEntries() < minc) fTip->Draw(); else fTip->Fit("gaus");
  fLip->Draw("same"); c1->cd();
  
  TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
  p5->SetFillColor(41); p5->SetFrameFillColor(10);
  fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
  TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
  TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
  hg->Divide(fFound,fGood,1.,1.,"B");
  hf->Divide(fFake,fGood,1.,1.,"B");
  hg->SetLineColor(4); hg->SetLineWidth(2);
  hg->SetMaximum(1.4);
  hg->SetYTitle("Tracking efficiency");
  hg->SetXTitle("Pt (GeV/c)");
  hg->Draw();
  
  TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
  line1->Draw("same");
  TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
  line2->Draw("same");
  
  hf->SetFillColor(1);
  hf->SetFillStyle(3013);
  hf->SetLineColor(2);
  hf->SetLineWidth(2);
  hf->Draw("histsame");
  TText* text = new TText(0.461176, 0.248448, "Fake tracks");
  text->SetTextSize(0.05);
  text->Draw();
  text = new TText(0.453919, 1.11408, "Good tracks");
  text->SetTextSize(0.05);
  text->Draw();
  
  TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
  
  TPad* p6 = new TPad("p6", "", 0., 0., 1., .5); p6->Draw();
  p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
  fE->SetFillColor(2); fE->SetFillStyle(3005);
  fE->SetXTitle("Arbitrary Units");
  if (fE->GetEntries() < minc) fE->Draw(); else fE->Fit("gaus"); c2->cd();
  
  TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
  p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
  fEp->SetFillColor(2); fEp->SetFillStyle(3005);
  fEp->SetMarkerStyle(8);
  fEp->SetMarkerSize(0.4);
  fEp->SetXTitle("p (GeV/c)"); fEp->SetYTitle("dE/dX (Arb. Units)");
  fEp->Draw(); c1->cd();
  
  TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
  c3->SetFillColor(42); c3->SetFrameFillColor(10);
  fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
  TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)", 
			 90, 0., 2.*TMath::Pi());
  hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B"); 
  hgphi->SetYTitle("Tracking efficiency");
  hgphi->SetXTitle("Phi (rad)");
  hgphi->Draw();
  
  TFile fc("AliITSComparison.root","RECREATE");
  c1->Write();
  c2->Write();
  c3->Write();
  fc.Close();
}
Commit	Line	Data
2d35f5d9	1	//-------------------------------------------------------------------------
	2	//
	3	// This is the PROOF-enabled version of ITS/AliITSComparisonV2.C macro.
	4	// Origin: Andrei.Zalite@cern.ch
	5	//
	6	//-------------------------------------------------------------------------
	7
f2247c19	8	#include "TChain.h"
f2247c19	9	#include "TH1F.h"
	10	#include "TH2F.h"
	11	#include "TList.h"
	12	#include "TClonesArray.h"
	13	#include "TCanvas.h"
	14	#include "TStyle.h"
	15	#include "TLine.h"
	16	#include "TText.h"
	17	#include "TFile.h"
f2247c19	18
	19	#include "AliLog.h"
	20	#include "AliVEvent.h"
	21	#include "AliESDEvent.h"
	22	#include "AliMCEvent.h"
	23	#include "AliESDtrack.h"
	24	#include "AliTrackReference.h"
	25	#include "AliMCComparisonTrack.h"
	26
	27	#include "AliITSComparisonTask.h"
	28
	29	ClassImp(AliITSComparisonTask)
	30
3d0bdd89	31	extern TStyle *gStyle;
3d0bdd89	32
f2247c19	33	AliITSComparisonTask::AliITSComparisonTask()
	34	: AliAnalysisTaskSE("AliITSComaprisonTask"),
	35	fListOfHistos(0),
3d0bdd89	36	fGood(0),
	37	fFound(0),
	38	fFake(0),
	39	fP(0),
	40	fL(0),
	41	fPt(0),
	42	fTip(0),
	43	fE(0),
	44	fEp(0),
	45	fGoodPhi(0),
	46	fFoundPhi(0),
	47	fLip(0)
f2247c19	48	{
	49	// Default constructor
	50	AliInfo("Default constructor AliITSComparisonTask");
	51	// Define input and output slots here
	52	// Input slot #0 works with a TChain
	53	DefineInput(0, TChain::Class());
	54	// Output slot #1 TList
	55	DefineOutput(1, TList::Class());
	56	}
	57
	58
	59
	60	AliITSComparisonTask::AliITSComparisonTask(const char* name)
	61	: AliAnalysisTaskSE(name),
	62	fListOfHistos(0),
3d0bdd89	63	fGood(0),
	64	fFound(0),
	65	fFake(0),
	66	fP(0),
	67	fL(0),
	68	fPt(0),
	69	fTip(0),
	70	fE(0),
	71	fEp(0),
	72	fGoodPhi(0),
	73	fFoundPhi(0),
	74	fLip(0)
f2247c19	75	{
	76	// Constructor
	77	AliInfo("Constructor AliITSComparisonTask");
	78	// Define input and output slots here
	79	// Input slot #0 works with a TChain
	80	DefineInput(0, TChain::Class());
	81	// Output slot #1 TList
	82	DefineOutput(1, TList::Class());
	83	}
	84
	85
	86
	87	void AliITSComparisonTask::UserCreateOutputObjects()
	88	{
	89	// Create histograms
	90	// Called once
	91	AliInfo("AliITSComparisonTask::UserCreateOutputObjects");
	92	// Create output container
	93	fListOfHistos = new TList();
	94
3d0bdd89	95	fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
	96	fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
	97	fFake = new TH1F("fFake", "Pt for fake tracks", 34, 0.2, 7.0);
	98	fP = new TH1F("fP", "PHI resolution", 50, -20., 20.);
	99	fL = new TH1F("fL", "LAMBDA resolution", 50, -20., 20.);
	100	fPt = new TH1F("fPt", "Relative Pt resolution", 30, -10., 10.);
	101	fTip = new TH1F("fTip", "Transverse impact parameter", 40, -800., 800.);
	102	fE = new TH1F("fE", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
	103	fEp = new TH2F("fEp", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
	104	fGoodPhi = new TH1F("fGoodPhi", "Phi for good tracks", 90, 0., 2.*TMath::Pi());
	105	fFoundPhi = new TH1F("fFoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
	106	fLip = new TH1F("fLip","Longitudinal impact parameter", 40, -800., 800.);
f2247c19	107
f2247c19	108
3d0bdd89	109	fListOfHistos->Add(fGood);
	110	fListOfHistos->Add(fFound);
	111	fListOfHistos->Add(fFake);
	112	fListOfHistos->Add(fP);
	113	fListOfHistos->Add(fL);
	114	fListOfHistos->Add(fPt);
	115	fListOfHistos->Add(fTip);
	116	fListOfHistos->Add(fE);
	117	fListOfHistos->Add(fEp);
	118	fListOfHistos->Add(fGoodPhi);
	119	fListOfHistos->Add(fFoundPhi);
	120	fListOfHistos->Add(fLip);
f2247c19	121	}
	122
	123
	124	void AliITSComparisonTask::UserExec(Option_t *)
	125	{
	126	// Main loop
	127	// Called for each event
	128
	129	// MC information
	130	AliMCEvent* mcEvent = MCEvent();
	131	if (!mcEvent) {
	132	Printf("ERROR: Could not retrieve MC event");
	133	return;
	134	}
	135
	136	Int_t nt = 0;
	137
	138	TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
	139
	140	// Loop over all MC tracks and select "good" tracks
	141	for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
84933864	142	AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
f2247c19	143	if (!track) {
	144	Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
	145	continue;
	146	}
	147
	148	// Track selection
	149	if (track->Pt() < 0.2) continue;
	150	if (track->Pt() > 7.0) continue;
	151	if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
	152
	153	Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
	154	if (TMath::Sqrt(vxvx+vyvy) > 3.5) continue;
	155	if (TMath::Abs(vz) > 50.) continue;
	156
	157	// Loop over Track References
3d0bdd89	158	Bool_t labelTPC = kFALSE;
f2247c19	159	AliTrackReference* trackRef = 0;
3d0bdd89	160	UInt_t iTSLayerMap = 0;
f2247c19	161	for (Int_t iTrackRef = 0; iTrackRef < track->GetNumberOfTrackReferences(); iTrackRef++) {
	162	trackRef = track->GetTrackReference(iTrackRef);
	163	if(trackRef) {
	164	Int_t detectorId = trackRef->DetectorId();
	165	if (detectorId == AliTrackReference::kITS) {
3d0bdd89	166	Float_t radius = trackRef->R();
	167	if (radius > 2.5 && radius < 5.5)
	168	iTSLayerMap \|= 1;
	169	else if (radius > 5.5 && radius < 8.5)
	170	iTSLayerMap \|= 1 << 1;
	171	else if (radius > 13.5 && radius < 16.5)
	172	iTSLayerMap \|= 1 << 2;
	173	else if (radius > 22. && radius < 26.)
	174	iTSLayerMap \|= 1 << 3;
	175	else if (radius > 36. && radius < 41.)
	176	iTSLayerMap \|= 1 << 4;
	177	else if (radius > 42. && radius < 46.)
	178	iTSLayerMap \|= 1 << 5;
f2247c19	179	else {
3d0bdd89	180	Printf("Wrong radius %f ", radius);
f2247c19	181	return;
	182	}
	183
	184	}
	185	if (detectorId == AliTrackReference::kTPC) {
3d0bdd89	186	labelTPC = kTRUE;
f2247c19	187	break;
	188	}
	189	}
	190	} // track references loop
	191
	192	// Skip tracks that passed not all ITS layers
3d0bdd89	193	if (iTSLayerMap != 0x3F) continue;
f2247c19	194
f2247c19	195	// "Good" tracks
3d0bdd89	196	if (labelTPC) {
f2247c19	197	AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
	198	ref->SetLabel(iTracks);
	199	TParticle* particle = track->Particle();
	200	Int_t pdg = particle->GetPdgCode();
	201	ref->SetPDG(pdg);
	202	ref->SetPz(track->Pz());
3d0bdd89	203	Float_t pt = track->Pt();
	204	ref->SetPt(pt);
	205	fGood->Fill(pt);
f2247c19	206	Float_t phig = track->Phi();
f2247c19	207	ref->SetPhi(phig);
3d0bdd89	208	fGoodPhi->Fill(phig);
f2247c19	209	nt++;
	210	}
	211	} // track loop
	212
	213
	214	// ESD information
	215	AliVEvent* event = InputEvent();
	216	if (!event) {
	217	Printf("ERROR: Could not retrieve event");
	218	return;
	219	}
	220
	221	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
	222
	223	Bool_t iFound;
	224	Int_t nfound = 0;
	225	Int_t nfake = 0;
	226	Int_t nlost = 0;
	227
3d0bdd89	228	Int_t mcGoods = refs->GetEntriesFast();
f2247c19	229
f2247c19	230	// Loop over all "good" MC tracks
3d0bdd89	231	for (Int_t k = 0; k < mcGoods; k++) {
f2247c19	232	AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
3d0bdd89	233	Int_t mcLabel = ref->GetLabel();
f2247c19	234	Float_t ptg = ref->GetPt();
3d0bdd89	235	Float_t phiG = ref->GetPhi();
f2247c19	236	iFound = kFALSE;
	237	for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
	238	AliESDtrack* track = esd->GetTrack(iTrack);
	239	if (!track) {
	240	Printf("ERROR: Could not receive track %d", iTrack);
	241	continue;
	242	}
	243
	244	// if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;
	245	if (! track->IsOn(AliESDtrack::kITSrefit)) continue;
	246
3d0bdd89	247	Int_t label = track->GetLabel();
f2247c19	248
3d0bdd89	249	if (mcLabel == TMath::Abs(label)) {
3d0bdd89	250	if (mcLabel == label) {
f2247c19	251	nfound++;
3d0bdd89	252	fFound->Fill(ptg);
3d0bdd89	253	fFoundPhi->Fill(phiG);
f2247c19	254	}
	255	else {
	256	nfake++;
3d0bdd89	257	fFake->Fill(ptg);
f2247c19	258	}
	259	iFound = kTRUE;
	260
	261	Float_t phi = track->Phi();
	262	Double_t pt = track->Pt();
	263
	264	Float_t lam = TMath::ATan2(track->Pz(),pt);
3d0bdd89	265	Float_t pt1 = 1/pt;
f2247c19	266
f2247c19	267	Float_t phig = ref->GetPhi();
3d0bdd89	268	fP->Fill((phi - phig)*1000.);
f2247c19	269
f2247c19	270	Float_t lamg = TMath::ATan2(ref->GetPz(),ptg);
3d0bdd89	271	fL->Fill((lam - lamg)*1000.);
f2247c19	272
3d0bdd89	273	fPt->Fill((pt1 - 1/ptg)/(1/ptg)*100.);
f2247c19	274
f2247c19	275	Float_t d,z; track->GetImpactParameters(d,z);
3d0bdd89	276	fTip->Fill(10000.*d);
3d0bdd89	277	fLip->Fill(10000.*z);
f2247c19	278
	279	Float_t mom = pt/TMath::Cos(lam);
	280	Float_t dedx = track->GetTPCsignal();
3d0bdd89	281	fEp->Fill(mom, dedx, 1.);
f2247c19	282
	283	Int_t pdg = ref->GetPDG();
	284	if (TMath::Abs(pdg )== 211) { // pions
	285	if (mom > 0.4 && mom < 0.5) {
3d0bdd89	286	fE->Fill(dedx, 1.);
f2247c19	287	}
	288	}
	289
	290	break;
	291	}
	292	}
	293	if (!iFound) {
	294	nlost++;
	295	}
	296	}
	297
	298	Printf(" Results: " );
	299	Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
	300
	301	refs->Clear();
	302
	303	// Post output data.
	304	PostData(1, fListOfHistos);
	305	}
	306
	307	//________________________________________________________________________
	308	void AliITSComparisonTask::Terminate(Option_t *)
	309	{
	310
	311	// Draw result to the screen
	312	// Called once at the end of the query
	313	fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
	314	if (!fListOfHistos) {
	315	Printf("ERROR: fListOfHistos not available");
	316	return;
	317	}
	318
3d0bdd89	319	fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
	320	fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
	321	fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
	322	fP = dynamic_cast<TH1F*>(fListOfHistos->At(3));
	323	fL = dynamic_cast<TH1F*>(fListOfHistos->At(4));
	324	fPt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
	325	fTip = dynamic_cast<TH1F*>(fListOfHistos->At(6));
	326	fE = dynamic_cast<TH1F*>(fListOfHistos->At(7));
	327	fEp = dynamic_cast<TH2F*>(fListOfHistos->At(8));
	328	fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
	329	fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
	330	fLip = dynamic_cast<TH1F*>(fListOfHistos->At(11));
f2247c19	331
	332	gStyle->SetOptStat(111110);
	333	gStyle->SetOptFit(1);
	334
	335	TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
	336
	337	Int_t minc = 33;
	338
	339	TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
	340	p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
3d0bdd89	341	fP->SetFillColor(4); fP->SetXTitle("(mrad)");
3d0bdd89	342	if (fP->GetEntries() < minc) fP->Draw(); else fP->Fit("gaus"); c1->cd();
f2247c19	343
	344	TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1., 0.6); p2->Draw();
	345	p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
3d0bdd89	346	fL->SetFillColor(4); fL->SetXTitle("(mrad)");
3d0bdd89	347	if (fL->GetEntries() < minc) fL->Draw(); else fL->Fit("gaus"); c1->cd();
f2247c19	348
	349	TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
	350	p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
3d0bdd89	351	fPt->SetFillColor(2); fPt->SetXTitle("(%)");
3d0bdd89	352	if (fPt->GetEntries() < minc) fPt->Draw(); else fPt->Fit("gaus"); c1->cd();
f2247c19	353
	354	TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
	355	p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
3d0bdd89	356	fTip->SetFillColor(6); fTip->SetXTitle("(micron)");
	357	if (fTip->GetEntries() < minc) fTip->Draw(); else fTip->Fit("gaus");
	358	fLip->Draw("same"); c1->cd();
f2247c19	359
	360	TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
	361	p5->SetFillColor(41); p5->SetFrameFillColor(10);
3d0bdd89	362	fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
f2247c19	363	TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
f2247c19	364	TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
3d0bdd89	365	hg->Divide(fFound,fGood,1.,1.,"B");
3d0bdd89	366	hf->Divide(fFake,fGood,1.,1.,"B");
f2247c19	367	hg->SetLineColor(4); hg->SetLineWidth(2);
	368	hg->SetMaximum(1.4);
	369	hg->SetYTitle("Tracking efficiency");
	370	hg->SetXTitle("Pt (GeV/c)");
	371	hg->Draw();
	372
	373	TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
	374	line1->Draw("same");
	375	TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
	376	line2->Draw("same");
	377
	378	hf->SetFillColor(1);
	379	hf->SetFillStyle(3013);
	380	hf->SetLineColor(2);
	381	hf->SetLineWidth(2);
	382	hf->Draw("histsame");
	383	TText* text = new TText(0.461176, 0.248448, "Fake tracks");
	384	text->SetTextSize(0.05);
	385	text->Draw();
	386	text = new TText(0.453919, 1.11408, "Good tracks");
	387	text->SetTextSize(0.05);
	388	text->Draw();
	389
	390	TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
	391
	392	TPad* p6 = new TPad("p6", "", 0., 0., 1., .5); p6->Draw();
	393	p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
3d0bdd89	394	fE->SetFillColor(2); fE->SetFillStyle(3005);
	395	fE->SetXTitle("Arbitrary Units");
	396	if (fE->GetEntries() < minc) fE->Draw(); else fE->Fit("gaus"); c2->cd();
f2247c19	397
	398	TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
	399	p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
3d0bdd89	400	fEp->SetFillColor(2); fEp->SetFillStyle(3005);
	401	fEp->SetMarkerStyle(8);
	402	fEp->SetMarkerSize(0.4);
	403	fEp->SetXTitle("p (GeV/c)"); fEp->SetYTitle("dE/dX (Arb. Units)");
	404	fEp->Draw(); c1->cd();
f2247c19	405
	406	TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
	407	c3->SetFillColor(42); c3->SetFrameFillColor(10);
3d0bdd89	408	fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
f2247c19	409	TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
f2247c19	410	90, 0., 2.*TMath::Pi());
3d0bdd89	411	hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B");
f2247c19	412	hgphi->SetYTitle("Tracking efficiency");
	413	hgphi->SetXTitle("Phi (rad)");
	414	hgphi->Draw();
	415
	416	TFile fc("AliITSComparison.root","RECREATE");
	417	c1->Write();
	418	c2->Write();
	419	c3->Write();
	420	fc.Close();
	421	}