Generation of D0 signal

[u/mrichter/AliRoot.git] / macros / CompareTime.C

//
// This macro draws time resolution for true hypothesis and fits it with gausian
// on a given reference plane in a given momentum range.
// The results is then printed to a PostScript file.
// The filename is composed out of input parameters.
//
// Input parameters:
// nRefPlane - indicate the reference plane.
//   1: in TPC  
//   2: out TPC
//   3: out TRD
// 
// minP, maxP - total momentum cuts. momentum from gAlice are used for cuts
// 
// The macro is using SORTED reference points in a file "trackRef.root".
// Track references can be sorted using CopyReferences.C
// 
// Before using check filenames and names of trees invlolved
// 
// ITS tracks (nRefPlane == 1) are compared at inner TPC ref Plane
// make shure they are properly propagated to this plane
//
// Sylwester Radomski, e-mail: S.Radomski@gsi.de
// Feb 14, 2002 
// 


CompareTime (Int_t nRefPlane, Double_t minP, Double_t maxP) 
{
  
  // Check input data Consistancy  
  if (nRefPlane < 1 || nRefPlane > 3) {
    cout << "Wrong Reference Plane Id ( " << nRefPalne << ") [1-3] " << endl;
    return;
  }

  if (maxP < minP) {
    cout << "MaxP lesser than MinP" << endl;
    return;
  }

  // Style

  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(1110);
  gStyle->SetOptFit(11);

  Double_t cc = 0.299792458;
  AliKalmanTrack::SetConvConst(100/cc/0.4);
  
  const char* gFileName = "galice.root";
  const char* refFileName = "trackRef.root";

  const char* trackFileName[] = {"AliTPCtracks.root", "AliTPCBackTracks.root", "AliTRDtracks.root"};
  const char *treeName[] = {"TreeT_ITSb_0", "seedsTPCtoTRD_0", "TRDb_0"};
  const char *refTreeName[] = {"TPC", "TPC", "TRD"};

  // Histograms

  const Int_t nTypes = 5; 
  Int_t codes[] = {11, 13, 211, 321, 2212};
  
  //Double_t cutoff[] = {20, 20, 20, 100, 100};
  Double_t cutoff[] = {30, 30, 30, 100, 200};

  const char *names[nTypes] = {"Electrons", "Muons", "Pions", "Kaons", "Protons"};
  
  TH1D *hDiffSame[nTypes];;

  for(Int_t i=0; i<nTypes; i++)
    hDiffSame[i] = new TH1D(names[i], names[i], 50, -cutoff[i], cutoff[i]);
  
  TH1D *hLength = new TH1D("length", "Length Resolution", 100, -0.5, 0.5);

  // Particles

  TFile *refFile = new TFile(gFileName, "READ");
  gAlice = (AliRun*)refFile->Get("gAlice");
  gAlice->GetEvent(0);


  // Reference tracks

  refFile = new TFile(refFileName, "READ");
  TTree *treeRef = (TTree*)refFile->Get("TreeTR0_Sorted"); 
  TBranch *trkRef = treeRef->GetBranch( refTreeName[nRefPlane-1] );
  TClonesArray *trkRefs = new TClonesArray("AliTrackReference", 100);
  trkRef->SetAddress(&trkRefs);
  
  Int_t ntracks = trkRef->GetEntries();
  cout << ntracks << endl;

  // Tracks
  
  TFile *trackFile = new TFile(trackFileName[nRefPlane-1]);

  TTree *tracks = (TTree*)trackFile->Get(treeName[nRefPlane-1]);
  TBranch *recTracks = tracks->GetBranch("tracks");
  AliTPCtrack *track = 0;
  recTracks->SetAddress(&track);

  for(Int_t t=0; t<tracks->GetEntries(); t++) {

    cout << t << "\r";
    
    recTracks->GetEvent(t);
    
    Int_t lab = track->GetLabel();
    if (lab < 0 || lab > 100000) continue;
    
    Int_t pdg = gAlice->Particle(lab)->GetPdgCode();
    treeRef->GetEvent(lab);
    
    Int_t nEntrRef = trkRefs->GetEntries();
    
    if (gAlice->Particle(lab)->P() < minP ) continue; 
    if (gAlice->Particle(lab)->P() > maxP ) continue;
    if (gAlice->Particle(lab)->GetFirstMother() != -1) continue;
    
    if (!(nEntrRef)) continue;
      
    Double_t timeTrue, timeTrack;
    Double_t trueLength;

    Int_t index =  nEntrRef-1;
    if (nRefPlane == 1) index = 0;

    timeTrue = ((AliTrackReference*)(*trkRefs)[index])->GetTime() * 1e12;
    trueLength = ((AliTrackReference*)(*trkRefs)[index])->GetLength();
    
    hLength->Fill(trueLength - track->GetIntegratedLength());
    
    for(Int_t i=0; i<nTypes; i++) {
      timeTrack = track->GetIntegratedTime(codes[i]);
      if (abs(pdg) == codes[i]) hDiffSame[i]->Fill(timeTrue - timeTrack);
    }
  }  
  
  TCanvas *c = new TCanvas();
  c->SetCanvasSize(640, 800);  
  c->SetWindowSize(650, 830);  
  c->Divide(2,3);
 
  c->cd(1);
  hLength->Draw();
  hLength->SetXTitle("True Length - Measured Length [cm]");
  
  for(Int_t i=0; i<nTypes; i++) {
    c->cd(2+i);
    gPad->SetGridx();
    gPad->SetGridy();
    hDiffSame[i]->Draw();  
    hDiffSame[i]->SetXTitle("time true - measured [ps]");    
    hDiffSame[i]->Fit("gaus");
  }

  char outFileName[80];
  const char *refPlaneNames[] = {"inTPC", "outTPC","outTRD"}; 
  
  c->Modified();
  c->ForceUpdate();

  sprintf(outFileName, "time_%s_%3.2f-%3.2f.ps", refPlaneNames[nRefPlane-1], minP, maxP);
  c->Print(outFileName);

  if (treeRef) delete treeRef;
  if (tracks) delete tracks;
  
}