[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDhistTRDpid.C

// ////////////////////////////////////////////////////////////////////////////
//   Macro to get TRD signal distribution from all six layers of TRD
//   This function takes the data from following directories 
//   "rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5",//   "rmom6", "rmom8", "rmom10"
//    corresponding to momenta
//   0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0
//   Running this macro will generate the root file TRDdEdxHistogramsV1.root 
//   containing histograms for all momenta.
//   Prashant Shukla (shukla@physi.uni-heidelberg.de)
// ////////////////////////////////////////////////////////////////////////////

#if !defined( __CINT__) || defined(__MAKECINT__)
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TPDGCode.h>
#include <TParticle.h>
#include <TTree.h>

#include "AliTRDtrack.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliESD.h"
#include "AliRun.h"
#include "AliStack.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#else
#endif


//_____________________________________________________________________________

Bool_t AliTRDhistTRDpid(Int_t oFile=1)
{
  //
  //

  const Int_t kNMom=11;

  TH1F *h1dEdxEL[kNMom];
  TH1F *h1dEdxPI[kNMom];
  TH1F *h1dEdxMU[kNMom];
  TH1F *h1dEdxKA[kNMom];
  TH1F *h1dEdxPR[kNMom];
  TH1F *h1MaxTimBinEL[kNMom];
  TH1F *h1MaxTimBinPI[kNMom];
  Char_t text[200];
 
  // Histograms for dedx of e and pi track in TRD 
  
  for (Int_t imom = 0; imom < kNMom; imom++) {
    sprintf(text,"h1dEdxEL%01d",imom+1);
    h1dEdxEL[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
    h1dEdxEL[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
    h1dEdxEL[imom]->GetYaxis()->SetTitle("Entries");
    h1dEdxEL[imom]->SetLineColor(kRed);
    
    sprintf(text,"h1dEdxPI%01d",imom+1);  
    h1dEdxPI[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
    //  TH1F *h1dEdxPI[imom] = new TH1F("","",250,0,2500);
    h1dEdxPI[imom]->GetXaxis()->SetTitle("dE/dx_{TRD} (arbitrary units)");
    h1dEdxPI[imom]->GetYaxis()->SetTitle("Entries");
    h1dEdxPI[imom]->SetLineColor(kBlue);

    sprintf(text,"h1dEdxMU%01d",imom+1);  
    h1dEdxMU[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
    h1dEdxMU[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
    h1dEdxMU[imom]->GetYaxis()->SetTitle("Entries");
    h1dEdxMU[imom]->SetLineColor(kGreen);
    
    sprintf(text,"h1dEdxKA%01d",imom+1);  
    h1dEdxKA[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
    h1dEdxKA[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
    h1dEdxKA[imom]->GetYaxis()->SetTitle("Entries");
    h1dEdxKA[imom]->SetLineColor(7);
    
    sprintf(text,"h1dEdxPR%01d",imom+1); 
    h1dEdxPR[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
    h1dEdxPR[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
    h1dEdxPR[imom]->GetYaxis()->SetTitle("Entries");
    h1dEdxPR[imom]->SetLineColor(6);
    
    sprintf(text,"h1MaxTimBinEL%01d",imom+1);  
    h1MaxTimBinEL[imom] = new TH1F(text,"Dist. Time Bin of max. Cluster for e(Red) and pi(Blue)",30,0,30);
    h1MaxTimBinEL[imom]->GetXaxis()->SetTitle("Time Bin of Maximum Cluster");
    h1MaxTimBinEL[imom]->GetYaxis()->SetTitle("Entries");
    h1MaxTimBinEL[imom]->SetLineColor(2);
    
    sprintf(text,"h1MaxTimBinPI%01d",imom+1);  
    h1MaxTimBinPI[imom] = new TH1F(text,"Time Bin of max. Cluster Pion",30,0,30);
    h1MaxTimBinPI[imom]->SetLineColor(4);
  }
  
  // PID
  Int_t partCode[5] = 
    {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
  
  // Files

  Char_t *FileDir[11] = {"rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5", "rmom6", "rmom8", "rmom10"}; 

  Float_t mome[11]={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0};

  Char_t gAliceFileName[200];
  Char_t esdFileName[200];
  for (Int_t imom = 0; imom < kNMom; imom++) {

    sprintf(gAliceFileName,"~/work/data/%s/galice.root", FileDir[imom]);
    sprintf(esdFileName,"~/work/data/%s/AliESDs.root", FileDir[imom]);

  // open run loader and load gAlice, kinematics and header
  AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName);
  if (!runLoader) {
    Error("CheckESD", "getting run loader from file %s failed",gAliceFileName);
    return kFALSE;
  }
  runLoader->LoadgAlice();
  gAlice = runLoader->GetAliRun();
  if (!gAlice) {
    Error("CheckESD", "no galice object found");
    return kFALSE;
  }
  runLoader->LoadKinematics();
  runLoader->LoadHeader();
  
  //  AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
  // open the ESD file
  TFile* esdFile = TFile::Open(esdFileName);
  if (!esdFile || !esdFile->IsOpen()) {
    Error("CheckESD", "opening ESD file %s failed",esdFileName);
    return kFALSE;
  }
  AliESD* esd = new AliESD;
  TTree* tree = (TTree*) esdFile->Get("esdTree");
  if (!tree) {
    Error("CheckESD", "no ESD tree found");
    return kFALSE;
  }
  tree->SetBranchAddress("ESD", &esd);

  // loop over events
  for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
    runLoader->GetEvent(iEvent);
    AliStack* stack = gAlice->Stack();
    Int_t nParticles = stack->GetNtrack();
    TArrayF vertex(3);
    runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
    Int_t nGene=0, nGenpi=0;
    for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
      TParticle* particle = stack->Particle(iParticle);
      if (!particle) continue;
      if (particle->Pt() < 0.001) continue;
      //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
      TVector3 dVertex(particle->Vx() - vertex[0], 
        	       particle->Vy() - vertex[1],
        	       particle->Vz() - vertex[2]);
      if (dVertex.Mag() > 0.0001) continue;
      if (TMath::Abs(particle->GetPdgCode()) == partCode[0]) nGene++;
      if (TMath::Abs(particle->GetPdgCode()) == partCode[2]) nGenpi++;
    }  // end loop over particles
    printf("Gen Prim. el= %d, Gen Prim. pi= %d in Event No.= %d in File = %s \n", nGene, nGenpi, iEvent, FileDir[imom]);
    
    // get the event summary data
    tree->GetEvent(iEvent);
    if (!esd) {
      Error("CheckESD", "no ESD object found for event %d", iEvent);
      return kFALSE;
    }

    ///////////////////////////////////////////
    //    AliTRDpidESD::MakePID(esd);
    //    AliESDpid::MakePID(esd);
    ///////////////////////////////////////////

    // Initializing number of electrons and pions
    Int_t nne=0, nnpi=0;    
    // loop over tracks
    Int_t nTracks = esd->GetNumberOfTracks();
    printf("Found %d tracks. \n",nTracks);

    for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
      AliESDtrack* track = esd->GetTrack(iTrack);
      // select tracks of selected particles
      Int_t label = TMath::Abs(track->GetLabel());
      if (label > stack->GetNtrack()) continue;     // background
      TParticle* particle = stack->Particle(label);
      //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
      TVector3 dVertex(particle->Vx() - vertex[0], 
        	       particle->Vy() - vertex[1],
        	       particle->Vz() - vertex[2]);
      if (dVertex.Mag() > 0.0001) continue;
      if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
      if (track->GetConstrainedChi2() > 1e9) continue;
      Double_t p[3];
      track->GetConstrainedPxPyPz(p);
      TVector3 pTrack(p);
      //      Double_t mom = track->GetP();

      // PID/////////////////////////////////////////////////////
      if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
      Int_t iGen = 5;
      for (Int_t i = 0; i < 5; i++) {
        if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;
      }

      // Filling muon kan and proton//////////////////////////////
      Int_t Itot=0;
      if(track->GetTRDsignal()) Itot=1;
      if(!Itot) continue;

      // dE/dx and TEST TRD
      Double_t TRDsignal = track->GetTRDsignal();
      //      if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal); 
      //      if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal); 
      //      if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal);
      //      if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal);
      //      if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal);

      //Filling elcron and pions////////////////////////////////
      const Int_t kNPlane = 6;
      // Track Selection 
      Double_t Ise=1.;
      for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
        Ise*=track->GetTRDsignals(iPlane);
      }
      if(Ise==0) continue;
      
      for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
        Double_t TRDsignal = track->GetTRDsignals(iPlane);
        if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal);
        if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal);
        if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal);
        if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal);
        if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal);

        if (iGen == 0) h1MaxTimBinEL[imom]->Fill(track->GetTRDTimBin(iPlane));
        if (iGen == 2) h1MaxTimBinPI[imom]->Fill(track->GetTRDTimBin(iPlane));
      }
      
      //      Int_t TRDncls = track->GetTRDncls();
      if (iGen == 0) nne++;
      if (iGen == 2) nnpi++;
    } // end of loop over tracks
    printf("Electron Tracks = %d, Pion Tracks =  %d \n", nne, nnpi);
  }// end of loop over events
  delete esd;
  esdFile->Close();
  delete esdFile;
  
  runLoader->UnloadHeader();
  runLoader->UnloadKinematics();
  delete runLoader;
  } // Loop over files
  
    
  Char_t fileprint[200];
  
  // draw the histograms if not in batch mode

  new TCanvas;
  gPad->SetBorderMode(0);
  gPad->SetFillColor(0);
  //  gStyle->SetStatColor(0);
  h1dEdxPI[0]->DrawCopy();
  h1dEdxEL[0]->DrawCopy("SAME");
  //  h1dEdxMU[0]->DrawCopy("SAME");
  //  h1dEdxKA[0]->DrawCopy("SAME");
  //  h1dEdxPR[0]->DrawCopy("SAME");

  new TCanvas;
  gPad->SetBorderMode(0);
  gPad->SetFillColor(0);
  h1MaxTimBinEL[0]->DrawCopy();
  h1MaxTimBinPI[0]->DrawCopy("SAME");


  // write the output histograms to a file
  
  Char_t fileresponse[200];
  sprintf(fileresponse,"TRDdEdxHistogramsV%d.root",oFile);
  
  TFile* outputFile = TFile::Open(fileresponse, "recreate");
  if (!outputFile || !outputFile->IsOpen()) {
    Error("CheckESD", "opening output file check.root failed");
    return kFALSE;
  }

  for (Int_t imom = 0; imom < kNMom; imom++) {
    h1dEdxPI[imom]->Write();
    h1dEdxEL[imom]->Write();
    h1dEdxMU[imom]->Write();
    h1dEdxKA[imom]->Write();
    h1dEdxPR[imom]->Write();
    h1MaxTimBinEL[imom]->Write();
    h1MaxTimBinPI[imom]->Write();
  }

  delete [] h1dEdxPI;  
  delete [] h1dEdxEL;
  delete [] h1dEdxMU;
  delete [] h1dEdxKA;
  delete [] h1dEdxPR;
  delete [] h1MaxTimBinEL;
  delete [] h1MaxTimBinPI;
  delete outputFile;
  
  // result of check
  Info("GetPIDsignals", "GetPIDsignals was successfull");
  return kTRUE;
  
} // end of main program
Commit	Line	Data
8470c96e	1	// ////////////////////////////////////////////////////////////////////////////
	2	// Macro to get TRD signal distribution from all six layers of TRD
	3	// This function takes the data from following directories
	4	// "rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5",// "rmom6", "rmom8", "rmom10"
	5	// corresponding to momenta
	6	// 0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0
	7	// Running this macro will generate the root file TRDdEdxHistogramsV1.root
	8	// containing histograms for all momenta.
	9	// Prashant Shukla (shukla@physi.uni-heidelberg.de)
	10	// ////////////////////////////////////////////////////////////////////////////
	11
	12	#if !defined( __CINT__) \|\| defined(__MAKECINT__)
	13	#include <TFile.h>
	14	#include <TH1.h>
	15	#include <TH2.h>
	16	#include <TCanvas.h>
	17	#include <TVector3.h>
	18	#include <TPDGCode.h>
	19	#include <TParticle.h>
	20	#include <TTree.h>
	21
	22	#include "AliTRDtrack.h"
	23	#include "AliRunLoader.h"
	24	#include "AliLoader.h"
	25	#include "AliESD.h"
	26	#include "AliRun.h"
	27	#include "AliStack.h"
	28	#include "AliHeader.h"
	29	#include "AliGenEventHeader.h"
	30	#else
	31	#endif
	32
	33
	34	//_____________________________________________________________________________
	35
	36	Bool_t AliTRDhistTRDpid(Int_t oFile=1)
	37	{
	38	//
	39	//
	40
	41	const Int_t kNMom=11;
	42
	43	TH1F *h1dEdxEL[kNMom];
	44	TH1F *h1dEdxPI[kNMom];
	45	TH1F *h1dEdxMU[kNMom];
	46	TH1F *h1dEdxKA[kNMom];
	47	TH1F *h1dEdxPR[kNMom];
	48	TH1F *h1MaxTimBinEL[kNMom];
	49	TH1F *h1MaxTimBinPI[kNMom];
	50	Char_t text[200];
	51
	52	// Histograms for dedx of e and pi track in TRD
	53
	54	for (Int_t imom = 0; imom < kNMom; imom++) {
	55	sprintf(text,"h1dEdxEL%01d",imom+1);
	56	h1dEdxEL[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000);
	57	h1dEdxEL[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
	58	h1dEdxEL[imom]->GetYaxis()->SetTitle("Entries");
	59	h1dEdxEL[imom]->SetLineColor(kRed);
	60
	61	sprintf(text,"h1dEdxPI%01d",imom+1);
	62	h1dEdxPI[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000);
	63	// TH1F *h1dEdxPI[imom] = new TH1F("","",250,0,2500);
	64	h1dEdxPI[imom]->GetXaxis()->SetTitle("dE/dx_{TRD} (arbitrary units)");
65	h1dEdxPI[imom]->GetYaxis()->SetTitle("Entries");
66	h1dEdxPI[imom]->SetLineColor(kBlue);
67
68	sprintf(text,"h1dEdxMU%01d",imom+1);
69	h1dEdxMU[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000);
70	h1dEdxMU[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
71	h1dEdxMU[imom]->GetYaxis()->SetTitle("Entries");
72	h1dEdxMU[imom]->SetLineColor(kGreen);
73
74	sprintf(text,"h1dEdxKA%01d",imom+1);
75	h1dEdxKA[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000);
76	h1dEdxKA[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
77	h1dEdxKA[imom]->GetYaxis()->SetTitle("Entries");
78	h1dEdxKA[imom]->SetLineColor(7);
79
80	sprintf(text,"h1dEdxPR%01d",imom+1);
81	h1dEdxPR[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000);
82	h1dEdxPR[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
83	h1dEdxPR[imom]->GetYaxis()->SetTitle("Entries");
84	h1dEdxPR[imom]->SetLineColor(6);
85
86	sprintf(text,"h1MaxTimBinEL%01d",imom+1);
87	h1MaxTimBinEL[imom] = new TH1F(text,"Dist. Time Bin of max. Cluster for e(Red) and pi(Blue)",30,0,30);
88	h1MaxTimBinEL[imom]->GetXaxis()->SetTitle("Time Bin of Maximum Cluster");
89	h1MaxTimBinEL[imom]->GetYaxis()->SetTitle("Entries");
90	h1MaxTimBinEL[imom]->SetLineColor(2);
91
92	sprintf(text,"h1MaxTimBinPI%01d",imom+1);
93	h1MaxTimBinPI[imom] = new TH1F(text,"Time Bin of max. Cluster Pion",30,0,30);
94	h1MaxTimBinPI[imom]->SetLineColor(4);
95	}
96
97	// PID
98	Int_t partCode[5] =
99	{kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
100
101	// Files
102
103	Char_t *FileDir[11] = {"rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5", "rmom6", "rmom8", "rmom10"};
104
105	Float_t mome[11]={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0};
106
107	Char_t gAliceFileName[200];
108	Char_t esdFileName[200];
109	for (Int_t imom = 0; imom < kNMom; imom++) {
110
111	sprintf(gAliceFileName,"~/work/data/%s/galice.root", FileDir[imom]);
112	sprintf(esdFileName,"~/work/data/%s/AliESDs.root", FileDir[imom]);
113
114	// open run loader and load gAlice, kinematics and header
115	AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName);
116	if (!runLoader) {
117	Error("CheckESD", "getting run loader from file %s failed",gAliceFileName);
118	return kFALSE;
119	}
120	runLoader->LoadgAlice();
121	gAlice = runLoader->GetAliRun();
122	if (!gAlice) {
123	Error("CheckESD", "no galice object found");
124	return kFALSE;
125	}
126	runLoader->LoadKinematics();
127	runLoader->LoadHeader();
128
129	// AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
130	// open the ESD file
131	TFile* esdFile = TFile::Open(esdFileName);
132	if (!esdFile \|\| !esdFile->IsOpen()) {
133	Error("CheckESD", "opening ESD file %s failed",esdFileName);
134	return kFALSE;
135	}
136	AliESD* esd = new AliESD;
137	TTree* tree = (TTree*) esdFile->Get("esdTree");
138	if (!tree) {
139	Error("CheckESD", "no ESD tree found");
140	return kFALSE;
141	}
142	tree->SetBranchAddress("ESD", &esd);
143
144	// loop over events
145	for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
146	runLoader->GetEvent(iEvent);
147	AliStack* stack = gAlice->Stack();
148	Int_t nParticles = stack->GetNtrack();
149	TArrayF vertex(3);
150	runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
151	Int_t nGene=0, nGenpi=0;
152	for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
153	TParticle* particle = stack->Particle(iParticle);
154	if (!particle) continue;
155	if (particle->Pt() < 0.001) continue;
156	// if (TMath::Abs(particle->Eta()) > 0.3) continue;
157	TVector3 dVertex(particle->Vx() - vertex[0],
158	particle->Vy() - vertex[1],
159	particle->Vz() - vertex[2]);
160	if (dVertex.Mag() > 0.0001) continue;
161	if (TMath::Abs(particle->GetPdgCode()) == partCode[0]) nGene++;
162	if (TMath::Abs(particle->GetPdgCode()) == partCode[2]) nGenpi++;
163	} // end loop over particles
164	printf("Gen Prim. el= %d, Gen Prim. pi= %d in Event No.= %d in File = %s \n", nGene, nGenpi, iEvent, FileDir[imom]);
165
166	// get the event summary data
167	tree->GetEvent(iEvent);
168	if (!esd) {
169	Error("CheckESD", "no ESD object found for event %d", iEvent);
170	return kFALSE;
171	}
172
173	///////////////////////////////////////////
174	// AliTRDpidESD::MakePID(esd);
175	// AliESDpid::MakePID(esd);
176	///////////////////////////////////////////
177
178	// Initializing number of electrons and pions
179	Int_t nne=0, nnpi=0;
180	// loop over tracks
181	Int_t nTracks = esd->GetNumberOfTracks();
182	printf("Found %d tracks. \n",nTracks);
183
184	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
185	AliESDtrack* track = esd->GetTrack(iTrack);
186	// select tracks of selected particles
187	Int_t label = TMath::Abs(track->GetLabel());
188	if (label > stack->GetNtrack()) continue; // background
189	TParticle* particle = stack->Particle(label);
190	// if (TMath::Abs(particle->Eta()) > 0.3) continue;
191	TVector3 dVertex(particle->Vx() - vertex[0],
192	particle->Vy() - vertex[1],
193	particle->Vz() - vertex[2]);
194	if (dVertex.Mag() > 0.0001) continue;
195	if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
196	if (track->GetConstrainedChi2() > 1e9) continue;
197	Double_t p[3];
198	track->GetConstrainedPxPyPz(p);
199	TVector3 pTrack(p);
200	// Double_t mom = track->GetP();
201
202	// PID/////////////////////////////////////////////////////
203	if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
204	Int_t iGen = 5;
205	for (Int_t i = 0; i < 5; i++) {
206	if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;
207	}
208
209	// Filling muon kan and proton//////////////////////////////
210	Int_t Itot=0;
211	if(track->GetTRDsignal()) Itot=1;
212	if(!Itot) continue;
213
214	// dE/dx and TEST TRD
215	Double_t TRDsignal = track->GetTRDsignal();
216	// if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal);
217	// if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal);
218	// if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal);
219	// if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal);
220	// if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal);
221
222	//Filling elcron and pions////////////////////////////////
223	const Int_t kNPlane = 6;
224	// Track Selection
225	Double_t Ise=1.;
226	for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
227	Ise*=track->GetTRDsignals(iPlane);
228	}
229	if(Ise==0) continue;
230
231	for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
232	Double_t TRDsignal = track->GetTRDsignals(iPlane);
233	if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal);
234	if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal);
235	if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal);
236	if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal);
237	if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal);
238
239	if (iGen == 0) h1MaxTimBinEL[imom]->Fill(track->GetTRDTimBin(iPlane));
240	if (iGen == 2) h1MaxTimBinPI[imom]->Fill(track->GetTRDTimBin(iPlane));
241	}
242
243	// Int_t TRDncls = track->GetTRDncls();
244	if (iGen == 0) nne++;
245	if (iGen == 2) nnpi++;
246	} // end of loop over tracks
247	printf("Electron Tracks = %d, Pion Tracks = %d \n", nne, nnpi);
248	}// end of loop over events
249	delete esd;
250	esdFile->Close();
251	delete esdFile;
252
253	runLoader->UnloadHeader();
254	runLoader->UnloadKinematics();
255	delete runLoader;
256	} // Loop over files
257
258
259	Char_t fileprint[200];
260
261	// draw the histograms if not in batch mode
262
263	new TCanvas;
264	gPad->SetBorderMode(0);
265	gPad->SetFillColor(0);
266	// gStyle->SetStatColor(0);
267	h1dEdxPI[0]->DrawCopy();
268	h1dEdxEL[0]->DrawCopy("SAME");
269	// h1dEdxMU[0]->DrawCopy("SAME");
270	// h1dEdxKA[0]->DrawCopy("SAME");
271	// h1dEdxPR[0]->DrawCopy("SAME");
272
273	new TCanvas;
274	gPad->SetBorderMode(0);
275	gPad->SetFillColor(0);
276	h1MaxTimBinEL[0]->DrawCopy();
277	h1MaxTimBinPI[0]->DrawCopy("SAME");
278
279
280	// write the output histograms to a file
281
282	Char_t fileresponse[200];
283	sprintf(fileresponse,"TRDdEdxHistogramsV%d.root",oFile);
284
285	TFile* outputFile = TFile::Open(fileresponse, "recreate");
286	if (!outputFile \|\| !outputFile->IsOpen()) {
287	Error("CheckESD", "opening output file check.root failed");
288	return kFALSE;
289	}
290
291	for (Int_t imom = 0; imom < kNMom; imom++) {
292	h1dEdxPI[imom]->Write();
293	h1dEdxEL[imom]->Write();
294	h1dEdxMU[imom]->Write();
295	h1dEdxKA[imom]->Write();
296	h1dEdxPR[imom]->Write();
297	h1MaxTimBinEL[imom]->Write();
298	h1MaxTimBinPI[imom]->Write();
299	}
300
301	delete [] h1dEdxPI;
302	delete [] h1dEdxEL;
303	delete [] h1dEdxMU;
304	delete [] h1dEdxKA;
305	delete [] h1dEdxPR;
306	delete [] h1MaxTimBinEL;
307	delete [] h1MaxTimBinPI;
308	delete outputFile;
309
310	// result of check
311	Info("GetPIDsignals", "GetPIDsignals was successfull");
312	return kTRUE;
313
314	} // end of main program
315