[u/mrichter/AliRoot.git] / ITS / ShowITSHitsRecPoints.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TCanvas.h>
#include <TClassTable.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TGeoManager.h>
#include <TH1.h>
#include <TInterpreter.h>
#include "AliCDBManager.h"
#include "AliGeomManager.h"
#include "AliHeader.h"
#include "AliITS.h"
#include "AliITSDetTypeRec.h"
#include "AliITSgeom.h"
#include "AliITSRecPoint.h"
#include "AliRun.h"
#endif

Int_t ShowITSHitsRecPoints(Bool_t align=kFALSE,
			   TString alignfile="ITSfullv11Misalignment.root")
{
  ///////////////////////////////////////////////////////////////////////
  // Macro to check clusters and hits in the 6 ITS layers              //
  ///////////////////////////////////////////////////////////////////////

  if (gClassTable->GetID("AliRun") < 0) {
    gInterpreter->ExecuteMacro("loadlibs.C");
  }
  else { 
    if(gAlice){
      delete AliRunLoader::Instance();
      delete gAlice;
      gAlice=0;
    }
  }
  // Set OCDB if needed
  AliCDBManager* man = AliCDBManager::Instance();
  if (!man->IsDefaultStorageSet()) {
    printf("Setting a local default storage and run number 0\n");
    man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
    man->SetRun(0);
  }
  else {
    printf("Using deafult storage \n");
  }
 
  // retrives geometry 
  if(!gGeoManager){
    AliGeomManager::LoadGeometry("geometry.root");
  }

  if(align) {
    TFile f(alignfile.Data());
    TClonesArray* ar = (TClonesArray*)f.Get("ITSAlignObjs");
    AliGeomManager::ApplyAlignObjsToGeom(*ar);
    f.Close();
  } else {
    AliGeomManager::ApplyAlignObjsFromCDB("ITS");
  }

  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  if (rl == 0x0){
    cerr<<"Can not open session RL=NULL"<< endl;
    return -1;
  }
  Int_t retval = rl->LoadgAlice();
  if (retval){
    cerr<<"LoadgAlice returned error"<<endl;
    return -1;
  }
  gAlice=rl->GetAliRun();

  retval = rl->LoadHeader();
  if (retval){
    cerr<<"LoadHeader returned error"<<endl;
    return -1;
  }


  AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  if(!ITSloader){
    cerr<<"ITS loader not found"<<endl;
    return -1;
  }
  ITSloader->LoadRecPoints("read");
  ITSloader->LoadHits("read");


  Float_t cluglo[3]={0.,0.,0.}; 
  AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
  ITS->SetTreeAddress();
  AliITSgeom *geom = ITS->GetITSgeom();
  AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
  detTypeRec->SetDefaults();

  Int_t modmin=geom->GetStartDet(0);
  Int_t modmax=geom->GetLastDet(2);
  Int_t totmod=modmax-modmin;
  Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  Float_t zlim[6]={15.,15.,22.,30.,45.,55.};

  TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
  TH1F** hmod=new TH1F*[6];
  TH1F** hxl=new TH1F*[6];
  TH1F** hzl=new TH1F*[6];
  TH1F** hxg=new TH1F*[6];
  TH1F** hyg=new TH1F*[6];
  TH2F** hxyg=new TH2F*[6];
  TH2F** hxygHits=new TH2F*[6];
  TH1F** hzg=new TH1F*[6];
  TH1F** hr=new TH1F*[6];
  TH1F** hphi=new TH1F*[6];
  TH1F** hq=new TH1F*[6];

  Char_t name[10];
  for(Int_t iLay=0;iLay<6;iLay++){
    sprintf(name,"hmod%d",iLay+1);
    hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
    hmod[iLay]->GetXaxis()->SetTitle("Module");
    hmod[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxloc%d",iLay+1);
    hxl[iLay]=new TH1F(name,"",100,-4.,4.);
    hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
    hxl[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hzloc%d",iLay+1);
    hzl[iLay]=new TH1F(name,"",100,-4.,4.);
    hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
    hzl[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxgl%d",iLay+1);
    hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
    hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
    hxg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hygl%d",iLay+1);
    hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
    hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
    hyg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxygl%d",iLay+1);
    hxyg[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
    hxyg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
    hxyg[iLay]->GetYaxis()->SetTitle("Yglob (cm)");
    sprintf(name,"hxygHitsl%d",iLay+1);
    hxygHits[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
    sprintf(name,"hzgl%d",iLay+1);
    hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
    hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
    hzg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hr%d",iLay+1);
    hr[iLay]=new TH1F(name,"",100,0.,50.);
    hr[iLay]->GetXaxis()->SetTitle("r (cm)");
    hr[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hphi%d",iLay+1);
    hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());    
    hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
    hphi[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hq%d",iLay+1);
    hq[iLay]=new TH1F(name,"",100,0.,300.);    
    hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
    hq[iLay]->GetXaxis()->CenterTitle();
  }

  TGraph **gpts=new TGraph*[3];
  TGraph2D **gpts3d=new TGraph2D*[3];
  TGraph **gRP=new TGraph*[6];
  TGraph **gHend=new TGraph*[6];
  TGraph **gHstart=new TGraph*[6];
  for(Int_t i=0;i<3;i++){
    gpts[i]=new TGraph(0);
    gpts3d[i]=new TGraph2D(0);
  }
  for(Int_t i=0;i<6;i++){
    gRP[i]=new TGraph(0);
    gRP[i]->GetXaxis()->SetTitle("global x [cm]");
    gRP[i]->GetYaxis()->SetTitle("global y [cm]");
    gHend[i]=new TGraph(0);
    gHstart[i]=new TGraph(0);
  }
  Int_t totev=rl->GetNumberOfEvents();
  printf("Total Number of events = %d\n",totev);

  Int_t iRP[6]={0},iH[6]={0};

  for(Int_t iev=0;iev<totev;iev++){
    rl->GetEvent(iev);
    TTree *TR = ITSloader->TreeR();
    TClonesArray *ITSrec  = detTypeRec->RecPoints();
    TBranch *branch = 0;
    if(TR && ITSrec){
      branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
      if(branch)branch->SetAddress(&ITSrec);
    }
    TTree *TH = ITSloader->TreeH();
    TClonesArray *hits=new TClonesArray("AliITShit",10000);
    TH->SetBranchAddress("ITS",&hits);

    Int_t nparticles = rl->GetHeader()->GetNtrack();
    cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;


    Int_t ipt=0;
    for(Int_t subd=0;subd<3;subd++){

      Int_t first = geom->GetStartDet(subd);
      Int_t last = geom->GetLastDet(subd);

      for (Int_t mod=first; mod<=last; mod++){
	detTypeRec->ResetRecPoints();
	branch->GetEvent(mod);
	Int_t nrecp = ITSrec->GetEntries();
	if(nrecp>0){
	  for(Int_t irec=0;irec<nrecp;irec++) {
	    AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
	    Int_t lay=recp->GetLayer();
	    hlayer->Fill(lay);
	    recp->GetGlobalXYZ(cluglo);
	    Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]); 
	    Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
	    if(iev<3){
 	      gpts[iev]->SetPoint(ipt,cluglo[0],cluglo[1]);
 	      gpts3d[iev]->SetPoint(ipt,cluglo[0],cluglo[1],cluglo[2]);
	      ipt++;
	    }
	    gRP[lay]->SetPoint(iRP[lay],cluglo[0],cluglo[1]); iRP[lay]++; 
	    hmod[lay]->Fill(mod);
	    hzl[lay]->Fill(recp->GetDetLocalZ());
	    hxl[lay]->Fill(recp->GetDetLocalX());
	    hzg[lay]->Fill(cluglo[2]);
	    hyg[lay]->Fill(cluglo[1]);
	    hxyg[lay]->Fill(cluglo[0],cluglo[1]);
	    hxg[lay]->Fill(cluglo[0]);
	    hr[lay]->Fill(rad);
	    hphi[lay]->Fill(phi);
	    hq[lay]->Fill(recp->GetQ());
	  }
	}
      }
    }


    cout<<" Now read hits "<<endl;

    Int_t nentrHits=(Int_t)TH->GetEntries();
    for (Int_t i=0; i<nentrHits; i++) {
      TH->GetEvent(i);
      Int_t nhit=hits->GetEntriesFast();
      for (Int_t ih=0; ih<nhit; ih++) {
	AliITShit *h=(AliITShit*)hits->UncheckedAt(ih);
	if(h->StatusExiting()) {
	  Double_t xl,yl,zl,tl,xl0,yl0,zl0,tl0;
	  h->GetPositionL(xl,yl,zl,tl);
	  h->GetPositionL0(xl0,yl0,zl0,tl0);
	  //if(TMath::Abs(yl-yl0)<0.0290) continue;
	  hxygHits[h->GetLayer()-1]->Fill(h->GetXG(),h->GetYG());
	  gHend[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],h->GetXG(),h->GetYG()); 
	  Double_t x0,y0,z0,t0;
	  h->GetPositionG0(x0,y0,z0,t0);
	  gHstart[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],x0,y0); 
	  iH[h->GetLayer()-1]++;
	  //printf("layer %d hit length xy %f hit length yl %f\n",h->GetLayer()-1,TMath::Sqrt((x0-h->GetXG())*(x0-h->GetXG())+(y0-h->GetYG())*(y0-h->GetYG())),TMath::Abs(yl-yl0));
	} 
      }
    }
  }


  TCanvas **c=new TCanvas*[6];
  Char_t ctit[30];
  for(Int_t iLay=0;iLay<6;iLay++){
    sprintf(name,"can%d",iLay+1);
    sprintf(ctit,"Layer %d",iLay+1);
    c[iLay]=new TCanvas(name,ctit,1200,900);
    c[iLay]->Divide(3,3,0.001,0.001);
    c[iLay]->cd(1);
    hmod[iLay]->Draw();
    c[iLay]->cd(2);
    hxl[iLay]->Draw();
    c[iLay]->cd(3);
    hzl[iLay]->Draw();
    c[iLay]->cd(4);
    hxg[iLay]->Draw();
    c[iLay]->cd(5);
    hyg[iLay]->Draw();
    c[iLay]->cd(6);
    hzg[iLay]->Draw();
    c[iLay]->cd(7);
    hr[iLay]->Draw();
    c[iLay]->cd(8);
    hphi[iLay]->Draw();    
    c[iLay]->cd(9);
    hxyg[iLay]->Draw();    
  }

  TCanvas *cev0;
  cev0=new TCanvas("cev0","Event 0",600,600);
  gpts[0]->SetMarkerStyle(7);
  gpts[0]->SetTitle(0);
  gpts[0]->Draw("AP");

  TCanvas *cev1;
  cev1=new TCanvas("cev1","Event 1",600,600);
  gpts[1]->SetMarkerStyle(7);
  gpts[1]->SetTitle(0);
  gpts[1]->Draw("AP");

  TCanvas *cev2;
  cev2=new TCanvas("cev2","Event 2",600,600);
  gpts[2]->SetMarkerStyle(7);
  gpts[2]->SetTitle(0);
  gpts[2]->Draw("AP");

  TCanvas *chr = new TCanvas("chr","chr");
  chr->Divide(3,2);
  for(Int_t i=0;i<6;i++) {
    chr->cd(i+1);
    gHend[i]->SetMarkerStyle(7);
    gHend[i]->SetMarkerColor(1);
    gHend[i]->Draw("A,P");
    gHstart[i]->SetMarkerStyle(7);
    gHstart[i]->SetMarkerColor(4);
    gHstart[i]->Draw("P");
    gRP[i]->SetMarkerStyle(7);
    gRP[i]->SetMarkerColor(2);
    gRP[i]->Draw("P");
  }


  return 0;
}


Int_t ShowITSHitsRecPointsNtuple(Bool_t align=kFALSE,
				 TString alignfile="ITSfullv11Misalignment.root")
{
  ///////////////////////////////////////////////////////////////////////
  // Macro to check clusters and hits in the 6 ITS layers              //
  // Creates also ntuple                                               //
  ///////////////////////////////////////////////////////////////////////

  if (gClassTable->GetID("AliRun") < 0) {
    gInterpreter->ExecuteMacro("loadlibs.C");
  }
  else { 
    if(gAlice){
      delete AliRunLoader::Instance();
      delete gAlice;
      gAlice=0;
    }
  }
  // Set OCDB if needed
  AliCDBManager* man = AliCDBManager::Instance();
  if (!man->IsDefaultStorageSet()) {
    printf("Setting a local default storage and run number 0\n");
    man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
    man->SetRun(0);
  }
  else {
    printf("Using deafult storage \n");
  }
 
  // retrives geometry 
  if(!gGeoManager){
    AliGeomManager::LoadGeometry("geometry.root");
  }
  if(align) {
    TFile f(alignfile.Data());
    TClonesArray* ar = (TClonesArray*)f.Get("ITSAlignObjs");
    AliGeomManager::ApplyAlignObjsToGeom(*ar);
    f.Close();
  } else {
    AliGeomManager::ApplyAlignObjsFromCDB("ITS");
  }


  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  if (rl == 0x0){
    cerr<<"Can not open session RL=NULL"<< endl;
    return -1;
  }
  Int_t retval = rl->LoadgAlice();
  if (retval){
    cerr<<"LoadgAlice returned error"<<endl;
    return -1;
  }
  gAlice=rl->GetAliRun();

  retval = rl->LoadHeader();
  if (retval){
    cerr<<"LoadHeader returned error"<<endl;
    return -1;
  }


  AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  if(!ITSloader){
    cerr<<"ITS loader not found"<<endl;
    return -1;
  }
  ITSloader->LoadRecPoints("read");
  ITSloader->LoadHits("read");


  Float_t cluglo[3]={0.,0.,0.}; 
  AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
  ITS->SetTreeAddress();
  AliITSgeom *geom = ITS->GetITSgeom();
  AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
  detTypeRec->SetDefaults();

  Int_t modmin=geom->GetStartDet(0);
  Int_t modmax=geom->GetLastDet(2);
  Int_t totmod=modmax-modmin;
  Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  Float_t zlim[6]={15.,15.,22.,30.,45.,55.};

  TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
  TH1F** hmod=new TH1F*[6];
  TH1F** hxl=new TH1F*[6];
  TH1F** hzl=new TH1F*[6];
  TH1F** hxg=new TH1F*[6];
  TH1F** hyg=new TH1F*[6];
  TH2F** hxyg=new TH2F*[6];
  TH2F** hxygHits=new TH2F*[6];
  TH1F** hzg=new TH1F*[6];
  TH1F** hr=new TH1F*[6];
  TH1F** hphi=new TH1F*[6];
  TH1F** hq=new TH1F*[6];

  Char_t name[10];
  for(Int_t iLay=0;iLay<6;iLay++){
    sprintf(name,"hmod%d",iLay+1);
    hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
    hmod[iLay]->GetXaxis()->SetTitle("Module");
    hmod[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxloc%d",iLay+1);
    hxl[iLay]=new TH1F(name,"",100,-4.,4.);
    hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
    hxl[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hzloc%d",iLay+1);
    hzl[iLay]=new TH1F(name,"",100,-4.,4.);
    hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
    hzl[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxgl%d",iLay+1);
    hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
    hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
    hxg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hygl%d",iLay+1);
    hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
    hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
    hyg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxygl%d",iLay+1);
    hxyg[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
    hxyg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
    hxyg[iLay]->GetYaxis()->SetTitle("Yglob (cm)");
    sprintf(name,"hxygHitsl%d",iLay+1);
    hxygHits[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
    sprintf(name,"hzgl%d",iLay+1);
    hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
    hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
    hzg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hr%d",iLay+1);
    hr[iLay]=new TH1F(name,"",100,0.,50.);
    hr[iLay]->GetXaxis()->SetTitle("r (cm)");
    hr[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hphi%d",iLay+1);
    hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());    
    hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
    hphi[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hq%d",iLay+1);
    hq[iLay]=new TH1F(name,"",100,0.,300.);    
    hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
    hq[iLay]->GetXaxis()->CenterTitle();
  }

  TGraph **gpts=new TGraph*[3];
  TGraph2D **gpts3d=new TGraph2D*[3];
  TGraph **gRP=new TGraph*[6];
  TGraph **gHend=new TGraph*[6];
  TGraph **gHstart=new TGraph*[6];
  for(Int_t i=0;i<3;i++){
    gpts[i]=new TGraph(0);
    gpts3d[i]=new TGraph2D(0);
  }
  for(Int_t i=0;i<6;i++){
    gRP[i]=new TGraph(0);
    gRP[i]->GetXaxis()->SetTitle("global x [cm]");
    gRP[i]->GetYaxis()->SetTitle("global y [cm]");
    gHend[i]=new TGraph(0);
    gHstart[i]=new TGraph(0);
  }
  Int_t totev=rl->GetNumberOfEvents();
  printf("Total Number of events = %d\n",totev);

  Int_t iRP[6]={0},iH[6]={0};

  TNtuple *nt = new TNtuple("nt","ntuple","lay:mod:deltaxl:deltax:deltay:deltaz:phi:z:xl:hlength");

  for(Int_t iev=0;iev<totev;iev++){
    rl->GetEvent(iev);
    TTree *TR = ITSloader->TreeR();
    TClonesArray *ITSrec  = detTypeRec->RecPoints();
    TBranch *branch = 0;
    if(TR && ITSrec){
      branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
      if(branch)branch->SetAddress(&ITSrec);
    }
    TTree *TH = ITSloader->TreeH();
    TClonesArray *hits=new TClonesArray("AliITShit",10000);
    TH->SetBranchAddress("ITS",&hits);

    Int_t nparticles = rl->GetHeader()->GetNtrack();
    cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;


    Int_t ipt=0;
    for(Int_t subd=0;subd<3;subd++){

      Int_t first = geom->GetStartDet(subd);
      Int_t last = geom->GetLastDet(subd);

      for (Int_t mod=first; mod<=last; mod++){
	detTypeRec->ResetRecPoints();
	branch->GetEvent(mod);
	Int_t nrecp = ITSrec->GetEntries();
	if(nrecp>0){
	  for(Int_t irec=0;irec<nrecp;irec++) {
	    AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
	    Int_t lay=recp->GetLayer();
	    hlayer->Fill(lay);
	    recp->GetGlobalXYZ(cluglo);
	    Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]); 
	    Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
	    if(cluglo[2]>0) {gRP[lay]->SetPoint(iRP[lay],cluglo[0],cluglo[1]); iRP[lay]++;} 
	    hmod[lay]->Fill(mod);
	    hzl[lay]->Fill(recp->GetDetLocalZ());
	    hxl[lay]->Fill(recp->GetDetLocalX());
	    hzg[lay]->Fill(cluglo[2]);
	    hyg[lay]->Fill(cluglo[1]);
	    hxyg[lay]->Fill(cluglo[0],cluglo[1]);
	    hxg[lay]->Fill(cluglo[0]);
	    hr[lay]->Fill(rad);
	    hphi[lay]->Fill(phi);
	    hq[lay]->Fill(recp->GetQ());

	    Double_t hlength=0.;
	    for(Int_t jhits=0;jhits<TH->GetEntries();jhits++) {
	      TH->GetEvent(jhits);
	      Int_t nhit=hits->GetEntriesFast();
	      for (Int_t ih=0; ih<nhit; ih++) {
		AliITShit *h=(AliITShit*)hits->UncheckedAt(ih);
		if(h->GetTrack()!=recp->GetLabel(0)) continue;
		if(h->GetLayer()-1!=lay) continue;
		if(h->GetModule()!=mod) continue;
		//hxygHits[h->GetLayer()-1]->Fill(h->GetXG(),h->GetYG());
		//gHend[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],h->GetXG(),h->GetYG()); 
		Double_t x0,y0,z0,t0;
		Double_t x,y,z,t;
		h->GetPositionG0(x0,y0,z0,t0);
		h->GetPositionG(x,y,z,t);
		//gHstart[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],x0,y0); 
		//iH[h->GetLayer()-1]++;
		Double_t xl,yl,zl,tl,xl0,yl0,zl0,tl0;
		h->GetPositionL(xl,yl,zl,tl);
		h->GetPositionL0(xl0,yl0,zl0,tl0);
		hlength += TMath::Abs(yl-yl0);
		if(!h->StatusExiting()) continue;
		nt->Fill(lay,
			 mod,
			 1.e4*(0.5*(xl+xl0)-recp->GetDetLocalX()),
			 1.e4*(0.5*(x+x0)-cluglo[0]),1.e4*(0.5*(y+y0)-cluglo[1]),
			 1.e4*(0.5*(z+z0)-cluglo[2]),
			 phi,
			 cluglo[2],
			 recp->GetDetLocalX(),
			 hlength);
		hlength = 0.;
	      }
	    }


	  }
	}
      }
    }


  }


  TCanvas **c=new TCanvas*[6];
  Char_t ctit[30];
  for(Int_t iLay=0;iLay<6;iLay++){
    sprintf(name,"can%d",iLay+1);
    sprintf(ctit,"Layer %d",iLay+1);
    c[iLay]=new TCanvas(name,ctit,1200,900);
    c[iLay]->Divide(3,3,0.001,0.001);
    c[iLay]->cd(1);
    hmod[iLay]->Draw();
    c[iLay]->cd(2);
    hxl[iLay]->Draw();
    c[iLay]->cd(3);
    hzl[iLay]->Draw();
    c[iLay]->cd(4);
    hxg[iLay]->Draw();
    c[iLay]->cd(5);
    hyg[iLay]->Draw();
    c[iLay]->cd(6);
    hzg[iLay]->Draw();
    c[iLay]->cd(7);
    hr[iLay]->Draw();
    c[iLay]->cd(8);
    hphi[iLay]->Draw();    
    c[iLay]->cd(9);
    hxyg[iLay]->Draw();    
  }

  TCanvas *cev0;
  cev0=new TCanvas("cev0","Event 0",600,600);
  gpts[0]->SetMarkerStyle(7);
  gpts[0]->SetTitle(0);
  gpts[0]->Draw("AP");

  TCanvas *cev1;
  cev1=new TCanvas("cev1","Event 1",600,600);
  gpts[1]->SetMarkerStyle(7);
  gpts[1]->SetTitle(0);
  gpts[1]->Draw("AP");

  TCanvas *cev2;
  cev2=new TCanvas("cev2","Event 2",600,600);
  gpts[2]->SetMarkerStyle(7);
  gpts[2]->SetTitle(0);
  gpts[2]->Draw("AP");

  TFile *out= new TFile("ntHitsRecPoints.root","recreate");
  nt->Write();
  out->Close();

  return 0;
}
Commit	Line	Data
3fd01708	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <TCanvas.h>
	3	#include <TClassTable.h>
	4	#include <TGraph.h>
	5	#include <TGraph2D.h>
	6	#include <TGeoManager.h>
	7	#include <TH1.h>
	8	#include <TInterpreter.h>
	9	#include "AliCDBManager.h"
	10	#include "AliGeomManager.h"
	11	#include "AliHeader.h"
	12	#include "AliITS.h"
	13	#include "AliITSDetTypeRec.h"
	14	#include "AliITSgeom.h"
	15	#include "AliITSRecPoint.h"
	16	#include "AliRun.h"
	17	#endif
	18
00fe805a	19	Int_t ShowITSHitsRecPoints(Bool_t align=kFALSE,
3fd01708	20	TString alignfile="ITSfullv11Misalignment.root")
	21	{
	22	///////////////////////////////////////////////////////////////////////
	23	// Macro to check clusters and hits in the 6 ITS layers //
	24	///////////////////////////////////////////////////////////////////////
	25
	26	if (gClassTable->GetID("AliRun") < 0) {
	27	gInterpreter->ExecuteMacro("loadlibs.C");
	28	}
	29	else {
	30	if(gAlice){
33c3c91a	31	delete AliRunLoader::Instance();
3fd01708	32	delete gAlice;
	33	gAlice=0;
	34	}
	35	}
	36	// Set OCDB if needed
	37	AliCDBManager* man = AliCDBManager::Instance();
	38	if (!man->IsDefaultStorageSet()) {
	39	printf("Setting a local default storage and run number 0\n");
162637e4	40	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
3fd01708	41	man->SetRun(0);
	42	}
	43	else {
	44	printf("Using deafult storage \n");
	45	}
	46
	47	// retrives geometry
	48	if(!gGeoManager){
	49	AliGeomManager::LoadGeometry("geometry.root");
	50	}
443c185f	51
3fd01708	52	if(align) {
	53	TFile f(alignfile.Data());
	54	TClonesArray* ar = (TClonesArray*)f.Get("ITSAlignObjs");
	55	AliGeomManager::ApplyAlignObjsToGeom(*ar);
	56	f.Close();
443c185f	57	} else {
443c185f	58	AliGeomManager::ApplyAlignObjsFromCDB("ITS");
3fd01708	59	}
	60
	61	AliRunLoader* rl = AliRunLoader::Open("galice.root");
	62	if (rl == 0x0){
	63	cerr<<"Can not open session RL=NULL"<< endl;
	64	return -1;
	65	}
	66	Int_t retval = rl->LoadgAlice();
	67	if (retval){
	68	cerr<<"LoadgAlice returned error"<<endl;
	69	return -1;
	70	}
	71	gAlice=rl->GetAliRun();
	72
	73	retval = rl->LoadHeader();
	74	if (retval){
	75	cerr<<"LoadHeader returned error"<<endl;
	76	return -1;
	77	}
	78
	79
	80	AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	81	if(!ITSloader){
	82	cerr<<"ITS loader not found"<<endl;
	83	return -1;
	84	}
	85	ITSloader->LoadRecPoints("read");
	86	ITSloader->LoadHits("read");
	87
	88
	89	Float_t cluglo[3]={0.,0.,0.};
	90	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
	91	ITS->SetTreeAddress();
	92	AliITSgeom *geom = ITS->GetITSgeom();
	93	AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
	94	detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
	95	detTypeRec->SetDefaults();
	96
	97	Int_t modmin=geom->GetStartDet(0);
	98	Int_t modmax=geom->GetLastDet(2);
	99	Int_t totmod=modmax-modmin;
	100	Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
	101	Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
	102
	103	TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
	104	TH1F** hmod=new TH1F*[6];
	105	TH1F** hxl=new TH1F*[6];
	106	TH1F** hzl=new TH1F*[6];
	107	TH1F** hxg=new TH1F*[6];
	108	TH1F** hyg=new TH1F*[6];
	109	TH2F** hxyg=new TH2F*[6];
	110	TH2F** hxygHits=new TH2F*[6];
	111	TH1F** hzg=new TH1F*[6];
	112	TH1F** hr=new TH1F*[6];
	113	TH1F** hphi=new TH1F*[6];
	114	TH1F** hq=new TH1F*[6];
	115
	116	Char_t name[10];
	117	for(Int_t iLay=0;iLay<6;iLay++){
	118	sprintf(name,"hmod%d",iLay+1);
	119	hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
	120	hmod[iLay]->GetXaxis()->SetTitle("Module");
	121	hmod[iLay]->GetXaxis()->CenterTitle();
	122	sprintf(name,"hxloc%d",iLay+1);
123	hxl[iLay]=new TH1F(name,"",100,-4.,4.);
124	hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
125	hxl[iLay]->GetXaxis()->CenterTitle();
126	sprintf(name,"hzloc%d",iLay+1);
127	hzl[iLay]=new TH1F(name,"",100,-4.,4.);
128	hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
129	hzl[iLay]->GetXaxis()->CenterTitle();
130	sprintf(name,"hxgl%d",iLay+1);
131	hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
132	hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
133	hxg[iLay]->GetXaxis()->CenterTitle();
134	sprintf(name,"hygl%d",iLay+1);
135	hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
136	hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
137	hyg[iLay]->GetXaxis()->CenterTitle();
138	sprintf(name,"hxygl%d",iLay+1);
139	hxyg[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
140	hxyg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
141	hxyg[iLay]->GetYaxis()->SetTitle("Yglob (cm)");
142	sprintf(name,"hxygHitsl%d",iLay+1);
143	hxygHits[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
144	sprintf(name,"hzgl%d",iLay+1);
145	hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
146	hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
147	hzg[iLay]->GetXaxis()->CenterTitle();
148	sprintf(name,"hr%d",iLay+1);
149	hr[iLay]=new TH1F(name,"",100,0.,50.);
150	hr[iLay]->GetXaxis()->SetTitle("r (cm)");
151	hr[iLay]->GetXaxis()->CenterTitle();
152	sprintf(name,"hphi%d",iLay+1);
153	hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
154	hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
155	hphi[iLay]->GetXaxis()->CenterTitle();
156	sprintf(name,"hq%d",iLay+1);
157	hq[iLay]=new TH1F(name,"",100,0.,300.);
158	hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
159	hq[iLay]->GetXaxis()->CenterTitle();
160	}
161
162	TGraph *gpts=new TGraph[3];
163	TGraph2D *gpts3d=new TGraph2D[3];
164	TGraph *gRP=new TGraph[6];
165	TGraph *gHend=new TGraph[6];
166	TGraph *gHstart=new TGraph[6];
167	for(Int_t i=0;i<3;i++){
168	gpts[i]=new TGraph(0);
169	gpts3d[i]=new TGraph2D(0);
170	}
171	for(Int_t i=0;i<6;i++){
172	gRP[i]=new TGraph(0);
173	gRP[i]->GetXaxis()->SetTitle("global x [cm]");
174	gRP[i]->GetYaxis()->SetTitle("global y [cm]");
175	gHend[i]=new TGraph(0);
176	gHstart[i]=new TGraph(0);
177	}
178	Int_t totev=rl->GetNumberOfEvents();
179	printf("Total Number of events = %d\n",totev);
180
181	Int_t iRP[6]={0},iH[6]={0};
182
183	for(Int_t iev=0;iev<totev;iev++){
184	rl->GetEvent(iev);
185	TTree *TR = ITSloader->TreeR();
186	TClonesArray *ITSrec = detTypeRec->RecPoints();
187	TBranch *branch = 0;
188	if(TR && ITSrec){
189	branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
190	if(branch)branch->SetAddress(&ITSrec);
191	}
192	TTree *TH = ITSloader->TreeH();
193	TClonesArray *hits=new TClonesArray("AliITShit",10000);
194	TH->SetBranchAddress("ITS",&hits);
195
196	Int_t nparticles = rl->GetHeader()->GetNtrack();
197	cout<<"Event #"<<iev<<" #Particles="<<nparticles<<endl;
198
199
200	Int_t ipt=0;
201	for(Int_t subd=0;subd<3;subd++){
202
203	Int_t first = geom->GetStartDet(subd);
204	Int_t last = geom->GetLastDet(subd);
205
206	for (Int_t mod=first; mod<=last; mod++){
207	detTypeRec->ResetRecPoints();
208	branch->GetEvent(mod);
209	Int_t nrecp = ITSrec->GetEntries();
210	if(nrecp>0){
211	for(Int_t irec=0;irec<nrecp;irec++) {
212	AliITSRecPoint recp = (AliITSRecPoint)ITSrec->At(irec);
213	Int_t lay=recp->GetLayer();
214	hlayer->Fill(lay);
215	recp->GetGlobalXYZ(cluglo);
216	Float_t rad=TMath::Sqrt(cluglo[0]cluglo[0]+cluglo[1]cluglo[1]);
217	Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
218	if(iev<3){
219	gpts[iev]->SetPoint(ipt,cluglo[0],cluglo[1]);
220	gpts3d[iev]->SetPoint(ipt,cluglo[0],cluglo[1],cluglo[2]);
221	ipt++;
222	}
223	gRP[lay]->SetPoint(iRP[lay],cluglo[0],cluglo[1]); iRP[lay]++;
224	hmod[lay]->Fill(mod);
225	hzl[lay]->Fill(recp->GetDetLocalZ());
226	hxl[lay]->Fill(recp->GetDetLocalX());
227	hzg[lay]->Fill(cluglo[2]);
228	hyg[lay]->Fill(cluglo[1]);
229	hxyg[lay]->Fill(cluglo[0],cluglo[1]);
230	hxg[lay]->Fill(cluglo[0]);
231	hr[lay]->Fill(rad);
232	hphi[lay]->Fill(phi);
233	hq[lay]->Fill(recp->GetQ());
234	}
235	}
236	}
237	}
238
239
240	cout<<" Now read hits "<<endl;
241
242	Int_t nentrHits=(Int_t)TH->GetEntries();
243	for (Int_t i=0; i<nentrHits; i++) {
244	TH->GetEvent(i);
245	Int_t nhit=hits->GetEntriesFast();
246	for (Int_t ih=0; ih<nhit; ih++) {
247	AliITShit h=(AliITShit)hits->UncheckedAt(ih);
248	if(h->StatusExiting()) {
249	Double_t xl,yl,zl,tl,xl0,yl0,zl0,tl0;
250	h->GetPositionL(xl,yl,zl,tl);
251	h->GetPositionL0(xl0,yl0,zl0,tl0);
252	//if(TMath::Abs(yl-yl0)<0.0290) continue;
253	hxygHits[h->GetLayer()-1]->Fill(h->GetXG(),h->GetYG());
254	gHend[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],h->GetXG(),h->GetYG());
255	Double_t x0,y0,z0,t0;
256	h->GetPositionG0(x0,y0,z0,t0);
257	gHstart[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],x0,y0);
258	iH[h->GetLayer()-1]++;
259	//printf("layer %d hit length xy %f hit length yl %f\n",h->GetLayer()-1,TMath::Sqrt((x0-h->GetXG())(x0-h->GetXG())+(y0-h->GetYG())(y0-h->GetYG())),TMath::Abs(yl-yl0));
260	}
261	}
262	}
263	}
264
265
266	TCanvas *c=new TCanvas[6];
267	Char_t ctit[30];
268	for(Int_t iLay=0;iLay<6;iLay++){
269	sprintf(name,"can%d",iLay+1);
270	sprintf(ctit,"Layer %d",iLay+1);
271	c[iLay]=new TCanvas(name,ctit,1200,900);
272	c[iLay]->Divide(3,3,0.001,0.001);
273	c[iLay]->cd(1);
274	hmod[iLay]->Draw();
275	c[iLay]->cd(2);
276	hxl[iLay]->Draw();
277	c[iLay]->cd(3);
278	hzl[iLay]->Draw();
279	c[iLay]->cd(4);
280	hxg[iLay]->Draw();
281	c[iLay]->cd(5);
282	hyg[iLay]->Draw();
283	c[iLay]->cd(6);
284	hzg[iLay]->Draw();
285	c[iLay]->cd(7);
286	hr[iLay]->Draw();
287	c[iLay]->cd(8);
288	hphi[iLay]->Draw();
289	c[iLay]->cd(9);
290	hxyg[iLay]->Draw();
291	}
292
293	TCanvas *cev0;
294	cev0=new TCanvas("cev0","Event 0",600,600);
295	gpts[0]->SetMarkerStyle(7);
296	gpts[0]->SetTitle(0);
297	gpts[0]->Draw("AP");
298
299	TCanvas *cev1;
300	cev1=new TCanvas("cev1","Event 1",600,600);
301	gpts[1]->SetMarkerStyle(7);
302	gpts[1]->SetTitle(0);
303	gpts[1]->Draw("AP");
304
305	TCanvas *cev2;
306	cev2=new TCanvas("cev2","Event 2",600,600);
307	gpts[2]->SetMarkerStyle(7);
308	gpts[2]->SetTitle(0);
309	gpts[2]->Draw("AP");
310
311	TCanvas *chr = new TCanvas("chr","chr");
312	chr->Divide(3,2);
313	for(Int_t i=0;i<6;i++) {
314	chr->cd(i+1);
315	gHend[i]->SetMarkerStyle(7);
316	gHend[i]->SetMarkerColor(1);
317	gHend[i]->Draw("A,P");
318	gHstart[i]->SetMarkerStyle(7);
319	gHstart[i]->SetMarkerColor(4);
320	gHstart[i]->Draw("P");
321	gRP[i]->SetMarkerStyle(7);
322	gRP[i]->SetMarkerColor(2);
323	gRP[i]->Draw("P");
324	}
325
326
327	return 0;
328	}
329
330
00fe805a	331	Int_t ShowITSHitsRecPointsNtuple(Bool_t align=kFALSE,
3fd01708	332	TString alignfile="ITSfullv11Misalignment.root")
	333	{
	334	///////////////////////////////////////////////////////////////////////
	335	// Macro to check clusters and hits in the 6 ITS layers //
	336	// Creates also ntuple //
	337	///////////////////////////////////////////////////////////////////////
	338
	339	if (gClassTable->GetID("AliRun") < 0) {
	340	gInterpreter->ExecuteMacro("loadlibs.C");
	341	}
	342	else {
	343	if(gAlice){
33c3c91a	344	delete AliRunLoader::Instance();
3fd01708	345	delete gAlice;
	346	gAlice=0;
	347	}
	348	}
	349	// Set OCDB if needed
	350	AliCDBManager* man = AliCDBManager::Instance();
	351	if (!man->IsDefaultStorageSet()) {
	352	printf("Setting a local default storage and run number 0\n");
162637e4	353	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
3fd01708	354	man->SetRun(0);
	355	}
	356	else {
	357	printf("Using deafult storage \n");
	358	}
	359
	360	// retrives geometry
	361	if(!gGeoManager){
	362	AliGeomManager::LoadGeometry("geometry.root");
	363	}
	364	if(align) {
	365	TFile f(alignfile.Data());
	366	TClonesArray* ar = (TClonesArray*)f.Get("ITSAlignObjs");
	367	AliGeomManager::ApplyAlignObjsToGeom(*ar);
	368	f.Close();
443c185f	369	} else {
443c185f	370	AliGeomManager::ApplyAlignObjsFromCDB("ITS");
3fd01708	371	}
3fd01708	372
443c185f	373
3fd01708	374	AliRunLoader* rl = AliRunLoader::Open("galice.root");
	375	if (rl == 0x0){
	376	cerr<<"Can not open session RL=NULL"<< endl;
	377	return -1;
	378	}
	379	Int_t retval = rl->LoadgAlice();
	380	if (retval){
	381	cerr<<"LoadgAlice returned error"<<endl;
	382	return -1;
	383	}
	384	gAlice=rl->GetAliRun();
	385
	386	retval = rl->LoadHeader();
	387	if (retval){
	388	cerr<<"LoadHeader returned error"<<endl;
	389	return -1;
	390	}
	391
	392
	393	AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	394	if(!ITSloader){
	395	cerr<<"ITS loader not found"<<endl;
	396	return -1;
	397	}
	398	ITSloader->LoadRecPoints("read");
	399	ITSloader->LoadHits("read");
	400
	401
	402	Float_t cluglo[3]={0.,0.,0.};
	403	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
	404	ITS->SetTreeAddress();
	405	AliITSgeom *geom = ITS->GetITSgeom();
	406	AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
	407	detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
	408	detTypeRec->SetDefaults();
	409
	410	Int_t modmin=geom->GetStartDet(0);
	411	Int_t modmax=geom->GetLastDet(2);
	412	Int_t totmod=modmax-modmin;
	413	Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
	414	Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
	415
	416	TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
	417	TH1F** hmod=new TH1F*[6];
	418	TH1F** hxl=new TH1F*[6];
	419	TH1F** hzl=new TH1F*[6];
	420	TH1F** hxg=new TH1F*[6];
	421	TH1F** hyg=new TH1F*[6];
	422	TH2F** hxyg=new TH2F*[6];
	423	TH2F** hxygHits=new TH2F*[6];
	424	TH1F** hzg=new TH1F*[6];
	425	TH1F** hr=new TH1F*[6];
	426	TH1F** hphi=new TH1F*[6];
	427	TH1F** hq=new TH1F*[6];
	428
	429	Char_t name[10];
	430	for(Int_t iLay=0;iLay<6;iLay++){
	431	sprintf(name,"hmod%d",iLay+1);
	432	hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
	433	hmod[iLay]->GetXaxis()->SetTitle("Module");
	434	hmod[iLay]->GetXaxis()->CenterTitle();
	435	sprintf(name,"hxloc%d",iLay+1);
	436	hxl[iLay]=new TH1F(name,"",100,-4.,4.);
	437	hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
438	hxl[iLay]->GetXaxis()->CenterTitle();
439	sprintf(name,"hzloc%d",iLay+1);
440	hzl[iLay]=new TH1F(name,"",100,-4.,4.);
441	hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
442	hzl[iLay]->GetXaxis()->CenterTitle();
443	sprintf(name,"hxgl%d",iLay+1);
444	hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
445	hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
446	hxg[iLay]->GetXaxis()->CenterTitle();
447	sprintf(name,"hygl%d",iLay+1);
448	hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
449	hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
450	hyg[iLay]->GetXaxis()->CenterTitle();
451	sprintf(name,"hxygl%d",iLay+1);
452	hxyg[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
453	hxyg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
454	hxyg[iLay]->GetYaxis()->SetTitle("Yglob (cm)");
455	sprintf(name,"hxygHitsl%d",iLay+1);
456	hxygHits[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
457	sprintf(name,"hzgl%d",iLay+1);
458	hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
459	hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
460	hzg[iLay]->GetXaxis()->CenterTitle();
461	sprintf(name,"hr%d",iLay+1);
462	hr[iLay]=new TH1F(name,"",100,0.,50.);
463	hr[iLay]->GetXaxis()->SetTitle("r (cm)");
464	hr[iLay]->GetXaxis()->CenterTitle();
465	sprintf(name,"hphi%d",iLay+1);
466	hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
467	hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
468	hphi[iLay]->GetXaxis()->CenterTitle();
469	sprintf(name,"hq%d",iLay+1);
470	hq[iLay]=new TH1F(name,"",100,0.,300.);
471	hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
472	hq[iLay]->GetXaxis()->CenterTitle();
473	}
474
475	TGraph *gpts=new TGraph[3];
476	TGraph2D *gpts3d=new TGraph2D[3];
477	TGraph *gRP=new TGraph[6];
478	TGraph *gHend=new TGraph[6];
479	TGraph *gHstart=new TGraph[6];
480	for(Int_t i=0;i<3;i++){
481	gpts[i]=new TGraph(0);
482	gpts3d[i]=new TGraph2D(0);
483	}
484	for(Int_t i=0;i<6;i++){
485	gRP[i]=new TGraph(0);
486	gRP[i]->GetXaxis()->SetTitle("global x [cm]");
487	gRP[i]->GetYaxis()->SetTitle("global y [cm]");
488	gHend[i]=new TGraph(0);
489	gHstart[i]=new TGraph(0);
490	}
491	Int_t totev=rl->GetNumberOfEvents();
492	printf("Total Number of events = %d\n",totev);
493
494	Int_t iRP[6]={0},iH[6]={0};
495
496	TNtuple *nt = new TNtuple("nt","ntuple","lay:mod:deltaxl:deltax:deltay:deltaz:phi:z:xl:hlength");
497
498	for(Int_t iev=0;iev<totev;iev++){
499	rl->GetEvent(iev);
500	TTree *TR = ITSloader->TreeR();
501	TClonesArray *ITSrec = detTypeRec->RecPoints();
502	TBranch *branch = 0;
503	if(TR && ITSrec){
504	branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
505	if(branch)branch->SetAddress(&ITSrec);
506	}
507	TTree *TH = ITSloader->TreeH();
508	TClonesArray *hits=new TClonesArray("AliITShit",10000);
509	TH->SetBranchAddress("ITS",&hits);
510
511	Int_t nparticles = rl->GetHeader()->GetNtrack();
512	cout<<"Event #"<<iev<<" #Particles="<<nparticles<<endl;
513
514
515	Int_t ipt=0;
516	for(Int_t subd=0;subd<3;subd++){
517
518	Int_t first = geom->GetStartDet(subd);
519	Int_t last = geom->GetLastDet(subd);
520
521	for (Int_t mod=first; mod<=last; mod++){
522	detTypeRec->ResetRecPoints();
523	branch->GetEvent(mod);
524	Int_t nrecp = ITSrec->GetEntries();
525	if(nrecp>0){
526	for(Int_t irec=0;irec<nrecp;irec++) {
527	AliITSRecPoint recp = (AliITSRecPoint)ITSrec->At(irec);
528	Int_t lay=recp->GetLayer();
529	hlayer->Fill(lay);
530	recp->GetGlobalXYZ(cluglo);
531	Float_t rad=TMath::Sqrt(cluglo[0]cluglo[0]+cluglo[1]cluglo[1]);
532	Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
533	if(cluglo[2]>0) {gRP[lay]->SetPoint(iRP[lay],cluglo[0],cluglo[1]); iRP[lay]++;}
534	hmod[lay]->Fill(mod);
535	hzl[lay]->Fill(recp->GetDetLocalZ());
536	hxl[lay]->Fill(recp->GetDetLocalX());
537	hzg[lay]->Fill(cluglo[2]);
538	hyg[lay]->Fill(cluglo[1]);
539	hxyg[lay]->Fill(cluglo[0],cluglo[1]);
540	hxg[lay]->Fill(cluglo[0]);
541	hr[lay]->Fill(rad);
542	hphi[lay]->Fill(phi);
543	hq[lay]->Fill(recp->GetQ());
544
545	Double_t hlength=0.;
546	for(Int_t jhits=0;jhits<TH->GetEntries();jhits++) {
547	TH->GetEvent(jhits);
548	Int_t nhit=hits->GetEntriesFast();
549	for (Int_t ih=0; ih<nhit; ih++) {
550	AliITShit h=(AliITShit)hits->UncheckedAt(ih);
551	if(h->GetTrack()!=recp->GetLabel(0)) continue;
552	if(h->GetLayer()-1!=lay) continue;
553	if(h->GetModule()!=mod) continue;
554	//hxygHits[h->GetLayer()-1]->Fill(h->GetXG(),h->GetYG());
555	//gHend[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],h->GetXG(),h->GetYG());
556	Double_t x0,y0,z0,t0;
557	Double_t x,y,z,t;
558	h->GetPositionG0(x0,y0,z0,t0);
559	h->GetPositionG(x,y,z,t);
560	//gHstart[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],x0,y0);
561	//iH[h->GetLayer()-1]++;
562	Double_t xl,yl,zl,tl,xl0,yl0,zl0,tl0;
563	h->GetPositionL(xl,yl,zl,tl);
564	h->GetPositionL0(xl0,yl0,zl0,tl0);
565	hlength += TMath::Abs(yl-yl0);
566	if(!h->StatusExiting()) continue;
567	nt->Fill(lay,
568	mod,
569	1.e4(0.5(xl+xl0)-recp->GetDetLocalX()),
570	1.e4(0.5(x+x0)-cluglo[0]),1.e4(0.5(y+y0)-cluglo[1]),
571	1.e4(0.5(z+z0)-cluglo[2]),
572	phi,
573	cluglo[2],
574	recp->GetDetLocalX(),
575	hlength);
576	hlength = 0.;
577	}
578	}
579
580
581	}
582	}
583	}
584	}
585
586
587	}
588
589
590	TCanvas *c=new TCanvas[6];
591	Char_t ctit[30];
592	for(Int_t iLay=0;iLay<6;iLay++){
593	sprintf(name,"can%d",iLay+1);
594	sprintf(ctit,"Layer %d",iLay+1);
595	c[iLay]=new TCanvas(name,ctit,1200,900);
596	c[iLay]->Divide(3,3,0.001,0.001);
597	c[iLay]->cd(1);
598	hmod[iLay]->Draw();
599	c[iLay]->cd(2);
600	hxl[iLay]->Draw();
601	c[iLay]->cd(3);
602	hzl[iLay]->Draw();
603	c[iLay]->cd(4);
604	hxg[iLay]->Draw();
605	c[iLay]->cd(5);
606	hyg[iLay]->Draw();
607	c[iLay]->cd(6);
608	hzg[iLay]->Draw();
609	c[iLay]->cd(7);
610	hr[iLay]->Draw();
611	c[iLay]->cd(8);
612	hphi[iLay]->Draw();
613	c[iLay]->cd(9);
614	hxyg[iLay]->Draw();
615	}
616
617	TCanvas *cev0;
618	cev0=new TCanvas("cev0","Event 0",600,600);
619	gpts[0]->SetMarkerStyle(7);
620	gpts[0]->SetTitle(0);
621	gpts[0]->Draw("AP");
622
623	TCanvas *cev1;
624	cev1=new TCanvas("cev1","Event 1",600,600);
625	gpts[1]->SetMarkerStyle(7);
626	gpts[1]->SetTitle(0);
627	gpts[1]->Draw("AP");
628
629	TCanvas *cev2;
630	cev2=new TCanvas("cev2","Event 2",600,600);
631	gpts[2]->SetMarkerStyle(7);
632	gpts[2]->SetTitle(0);
633	gpts[2]->Draw("AP");
634
635	TFile *out= new TFile("ntHitsRecPoints.root","recreate");
636	nt->Write();
637	out->Close();
638
639	return 0;
640	}