[u/mrichter/AliRoot.git] / ITS / UPGRADE / testITSU / compClusHits.C

TObjArray histoArr;
enum {kNPixAll=0,kNPixSPL=1,kDR=0,kDTXodd,kDTXeven,kDTZ, kDTXoddSPL,kDTXevenSPL,kDTZSPL};

TPaveStats* GetStPad(TH1* hst);
TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl=-1,Int_t col=-1);
TCanvas* DrawNP(int np, TObjArray* harr=0, TCanvas* cnv=0);
TH1* GetHistoClSize(int npix,int id,TObjArray* harr=0);
void DrawReport(const char* psname, TObjArray* harr=0);


typedef struct {
  Int_t evID;
  Int_t volID;
  Int_t lrID;
  Int_t clID;
  Int_t nPix;
  Int_t nX;
  Int_t nZ;
  Float_t pt;
  Float_t eta;
  Float_t xyz[3];
  Float_t dX;
  Float_t dZ;  
  Bool_t split;  
  Bool_t prim;
  Int_t  pdg;
  Int_t  ntr;
} clSumm;

void compClusHits(int nev=-1)
{
  const int kSplit=0x1<<22;
  const int kSplCheck=0x1<<23;
  //
  gSystem->Load("libITSUpgradeBase");
  gSystem->Load("libITSUpgradeSim");
  gSystem->Load("libITSUpgradeRec");
  gROOT->SetStyle("Plain");

  gAlice=NULL;
  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  runLoader->LoadgAlice();

  gAlice = runLoader->GetAliRun();

  runLoader->LoadHeader();
  runLoader->LoadKinematics();
  runLoader->LoadRecPoints();
  runLoader->LoadSDigits();
  runLoader->LoadHits();

  AliLoader *dl = runLoader->GetDetectorLoader("ITS");

  AliGeomManager::LoadGeometry("geometry.root");
  AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
  Int_t nLayers = gm->GetNLayers();
  AliITSUClusterPix::SetGeom(gm);
  //
  AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
  its->CreateClusterArrays();
  //
  Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
  //
  TTree * cluTree = 0x0;
  TTree *hitTree = 0x0;
  TClonesArray *hitList=new TClonesArray("AliITSUHit");
  //
  TObjArray arrMCTracks; // array of hit arrays for each particle
  //
  Float_t xyzClGloF[3];
  Double_t xyzClGlo[3],xyzClTr[3];
  Int_t labels[3];
  int nLab = 0;
  int nlr=its->GetNLayersActive();
  int ntotev = (Int_t)runLoader->GetNumberOfEvents();
  printf("N Events : %i \n",ntotev);
  if (nev>0) ntotev = TMath::Min(nev,ntotev);
  //
  // output tree
  TFile* flOut = TFile::Open("clInfo.root","recreate");
  TTree* trOut = new TTree("clitsu","clitsu");
  clSumm cSum;
  trOut->Branch("evID", &cSum.evID ,"evID/I");
  trOut->Branch("volID",&cSum.volID,"volID/I");
  trOut->Branch("lrID", &cSum.lrID ,"lrID/I");  
  trOut->Branch("clID", &cSum.clID ,"clID/I");  
  trOut->Branch("nPix", &cSum.nPix ,"nPix/I");
  trOut->Branch("nX"  , &cSum.nX   ,"nX/I");
  trOut->Branch("nZ"  , &cSum.nZ   ,"nZ/I");
  trOut->Branch("pt"  , &cSum.pt   ,"pt/F");  
  trOut->Branch("eta"  ,&cSum.eta  ,"eta/F");  
  trOut->Branch("xyz",   cSum.xyz,  "xyz[3]/F");  
  trOut->Branch("dX"  , &cSum.dX   ,"dX/F");
  trOut->Branch("dZ"  , &cSum.dZ   ,"dZ/F");  
  trOut->Branch("split",&cSum.split,"split/O");
  trOut->Branch("prim", &cSum.prim, "prim/O");
  trOut->Branch("pdg",  &cSum.pdg,  "pdg/I");
  trOut->Branch("ntr",  &cSum.ntr,  "ntr/I");
  //
  for (Int_t iEvent = 0; iEvent < ntotev; iEvent++) {
    printf("\n Event %i \n",iEvent);
    runLoader->GetEvent(iEvent);
    AliStack *stack = runLoader->Stack();
    cluTree=dl->TreeR();
    hitTree=dl->TreeH();
    hitTree->SetBranchAddress("ITS",&hitList);
    // 
    // read clusters
    for (int ilr=nlr;ilr--;) {
      TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr));
      if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);}
      br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress());
    }
    cluTree->GetEntry(0); 
    its->ProcessClusters();
    //
    // read hits
    for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
      hitTree->GetEntry(iEnt);
      int nh = hitList->GetEntries();      
      for(Int_t iHit=0; iHit<hitList->GetEntries();iHit++){
	AliITSUHit *pHit = (AliITSUHit*)hitList->At(iHit);
	int mcID = pHit->GetTrack();
	TClonesArray* harr = arrMCTracks.GetEntriesFast()>mcID ? (TClonesArray*)arrMCTracks.At(mcID) : 0;
	if (!harr) {
	  harr = new TClonesArray("AliITSUHit"); // 1st encounter of the MC track
	  arrMCTracks.AddAtAndExpand(harr,mcID);
	}
	//
	new ( (*harr)[harr->GetEntriesFast()] ) AliITSUHit(*pHit);
      }
    }
    //
    // compare clusters and hits
    //
    printf(" tree entries: %d\n",cluTree->GetEntries());
    //
    for (int ilr=0;ilr<nlr;ilr++) {
      AliITSURecoLayer* lr = its->GetLayerActive(ilr);
      TClonesArray* clr = lr->GetClusters();
      int nClu = clr->GetEntries();
      printf("Layer %d : %d clusters\n",ilr,nClu);
      //
      for (int icl=0;icl<nClu;icl++) {
	AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(icl);
	int modID = cl->GetVolumeId();
	
	//------------ check if this is a split cluster
	int sInL = modID - gm->GetFirstModIndex(ilr);
	if (!cl->TestBit(kSplCheck)) {
	  cl->SetBit(kSplCheck);
	  // check if there is no other cluster with same label on this module
	  AliITSURecoSens* sens = lr->GetSensor(sInL);
	  int nclSn = sens->GetNClusters();
	  int offs = sens->GetFirstClusterId();
	  //	printf("To check for %d (mod:%d) N=%d from %d\n",icl,modID,nclSn,offs);
	  for (int ics=0;ics<nclSn;ics++) {
	    AliITSUClusterPix* clusT = (AliITSUClusterPix*)lr->GetCluster(offs+ics); // access to clusters
	    if (clusT==cl) continue;
	    for (int ilb0=0;ilb0<3;ilb0++) {
	      int lb0 = cl->GetLabel(ilb0); if (lb0<=-1) break;
	      for (int ilb1=0;ilb1<3;ilb1++) {
		int lb1 = clusT->GetLabel(ilb1); if (lb1<=-1) break;
		if (lb1==lb0) {
		  cl->SetBit(kSplit);
		  clusT->SetBit(kSplit);
		  /*
		    printf("Discard clusters of module %d:\n",modID);
		    cl->Print();
		    clusT->Print();
		  */
		  break;
		}
	      }
	    }
	  }
	}
	//------------
	AliITSsegmentation* segm = gm->GetSegmentation(ilr);
	//
	cl->GetGlobalXYZ(xyzClGloF);
	int clsize = cl->GetNPix();
	for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
	TGeoHMatrix* mat = gm->GetMatrixSens(modID);
	if (!mat) {printf("failed to get matrix for module %d\n",cl->GetVolumeId());}
	mat->MasterToLocal(xyzClGlo,xyzClTr);
	//
	int col,row;
	segm->LocalToDet(xyzClTr[0],xyzClTr[2],row,col); // effective col/row
	nLab = 0;
	for (int il=0;il<3;il++) {
	  if (cl->GetLabel(il)>=0) labels[nLab++] = cl->GetLabel(il);
	  else break;
	}
	// find hit info
	for (int il=0;il<nLab;il++) {
	  TClonesArray* htArr = (TClonesArray*)arrMCTracks.At(labels[il]);
	  if (!htArr) {printf("did not find MChits for label %d ",labels[il]); cl->Print(); continue;}
	  //
	  int nh = htArr->GetEntriesFast();
	  AliITSUHit *pHit=0;
	  for (int ih=nh;ih--;) {
	    AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih);
	    if (tHit->GetModule()!=modID) continue;
	    pHit = tHit;
	    break;
	  }
	  if (!pHit) {
	    printf("did not find MChit for label %d on module %d ",il,modID); 
	    cl->Print(); 
	    htAtt->Print();
	    continue;
	  }
	  //
	  pHit->GetPositionG(xg1,yg1,zg1);
	  pHit->GetPositionG0(xg0,yg0,zg0,tg0);
	  //
	  double txyzH[3],gxyzH[3] = { (xg1+xg0)/2, (yg1+yg0)/2, (zg1+zg0)/2 };
	  mat->MasterToLocal(gxyzH,txyzH);
	  double rcl = TMath::Sqrt(xyzClTr[0]*xyzClTr[0]+xyzClTr[1]*xyzClTr[1]);
	  double rht = TMath::Sqrt(txyzH[0]*txyzH[0]+txyzH[1]*txyzH[1]);
	  //
	  GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
	  if (cl->TestBit(kSplit)) {
	    if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
	    else       GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
	    GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
	    GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize);
	  }
	  if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
	  else       GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
	  GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
	  GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize);
	  //
	  cSum.evID = iEvent;
	  cSum.volID = cl->GetVolumeId();
	  cSum.lrID = ilr;
	  cSum.clID = icl;
	  cSum.nPix = cl->GetNPix();
	  cSum.nX   = cl->GetNx();
	  cSum.nZ   = cl->GetNz();
	  cSum.split = cl->TestBit(kSplit);
	  cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
	  cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
	  int label = cl->GetLabel(0);
	  TParticle* part = 0;
	  if (label>=0 && (part=stack->Particle(label)) ) {
	    cSum.pdg = part->GetPdgCode();
	    cSum.eta = part->Eta();
	    cSum.pt  = part->Pt();
	    cSum.prim = stack->IsPhysicalPrimary(label);
	  }
	  cSum.ntr = 0;
	  for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
	  for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
	  //
	  trOut->Fill();
	  /*
	  if (clsize==5) {
	    printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
		   modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
		   gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
	    cl->Print();
	    pHit->Print();
	    //
	    double a0,b0,c0,a1,b1,c1,e0;
	    pHit->GetPositionL0(a0,b0,c0,e0);
	    pHit->GetPositionL(a1,b1,c1);
	    float cloc[3];
	    cl->GetLocalXYZ(cloc);
	    printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
	    printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
	  }
	  */
	  //
	}
      }
    }
    
    //    layerClus.Clear();
    //
    arrMCTracks.Delete();
  }//event loop
  //
  flOut->cd();
  trOut->Write();
  delete trOut;
  flOut->Close();
  flOut->Delete();
  DrawReport("clinfo.ps",&histoArr);
  //
}

void DrawReport(const char* psname, TObjArray* harr) 
{
  gStyle->SetOptFit(1);
  if (!harr) harr = &histoArr;
  TCanvas* cnv = new TCanvas("cl","cl",900,600);
  //
  TString psnm1 = psname;
  if (psnm1.IsNull()) psnm1 = "clusters.ps";
  TString psnm0 = psnm1.Data(); 
  psnm0 += "[";
  TString psnm2 = psnm1.Data(); 
  psnm2 += "]";
  cnv->Print(psnm0.Data());
  //
  TH1* clall = GetHistoClSize(0,kNPixAll,harr);
  clall->SetLineColor(kRed);
  clall->Draw();
  TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr);
  clSpl->SetLineColor(kBlue);
  clSpl->Draw("sames");
  gPad->Modified();
  gPad->Update();
  SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor());
  SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor());
  gPad->Modified();
  gPad->Update();
  gPad->SetLogy(1);
  //
  cnv->cd();
  cnv->Print(psnm1.Data());
  //
  // plot cluster sized from 1 to 10
  for (int i=1;i<=10;i++) {
    if (clall->GetBinContent(clall->FindBin(i))<100) continue;
    DrawNP(i,harr,cnv);
    cnv->Print(psnm1.Data());
  }
  cnv->Print(psnm2.Data());
}

//------------------------------------------
TH1* GetHistoClSize(int npix,int id,TObjArray* harr)
{
  // book histos
  TH1* h = 0;
  if (!harr) harr = &histoArr;
  //
  if (npix<1) {
    if (harr->GetEntriesFast()>=id && (h=(TH1*)harr->At(id))) return h;
    h = new TH1F("npixAll","npixAll",150,0.5,54.5);
    h->SetLineColor(kRed);
    harr->AddAtAndExpand(h, kNPixAll);
    //
    h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
    h->SetLineColor(kBlue);
    harr->AddAtAndExpand(h, kNPixSPL);
    //
    h = (TH1*)harr->At(id);
    if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
    return h;
  }
  //
  int idh = npix*10+id;
  if (harr->GetEntriesFast()>=idh && (h=(TH1*)harr->At(idh))) return h;
  //
  const int nbin=100;
  const double kdiff=80;
  // need to create set of histos
  //
  h = new TH1F(Form("dxy_npix%d",npix),Form("dr_npix%d",npix),nbin,-kdiff,kdiff);
  harr->AddAtAndExpand(h, npix*10 + kDR);
  //
  h  = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
  h->SetLineColor(kRed);
  harr->AddAtAndExpand(h, npix*10 + kDTXodd);
  h  = new TH1F(Form("dtxEVN_npix%d",npix),Form("dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
  h->SetLineColor(kBlue);
  harr->AddAtAndExpand(h, npix*10 + kDTXeven);
  //
  h  = new TH1F(Form("dtz_npix%d",npix),Form("dtz_npix%d",npix),nbin,-kdiff,kdiff);
  h->SetLineColor(kGreen);
  harr->AddAtAndExpand(h, npix*10 + kDTZ);
  //
  //
  h  = new TH1F(Form("SPL_dtxODD_npix%d",npix),Form("SPL_dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
  h->SetLineColor(kMagenta);
  h->SetFillColor(kMagenta);
  h->SetFillStyle(3001);  
  h->SetLineStyle(2);
  harr->AddAtAndExpand(h, npix*10 + kDTXoddSPL);
  h  = new TH1F(Form("SPL_dtxEVN_npix%d",npix),Form("SPL_dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
  h->SetLineColor(kCyan);
  h->SetFillColor(kCyan);
  h->SetFillStyle(3006);  
  h->SetLineStyle(2);
  harr->AddAtAndExpand(h, npix*10 + kDTXevenSPL);
  //
  h  = new TH1F(Form("SPL_dtz_npix%d",npix),Form("SPLdtz_npix%d",npix),nbin,-kdiff,kdiff);
  harr->AddAtAndExpand(h, npix*10 + kDTZSPL);
  //
  h->SetLineColor(kGreen+2);
  h->SetFillColor(kGreen+2);
  h->SetLineStyle(2);
  h->SetFillStyle(3001);
  h = (TH1*)harr->At(idh);
  if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
  return h;
}


TCanvas* DrawNP(int np, TObjArray* harr, TCanvas* cnv)
{
  if (!harr) harr = &histoArr;
  if (!cnv) cnv = new TCanvas(Form("cnv%d",np),Form("cnv%d",np),900,700);
  cnv->Clear();
  cnv->Divide(2,1);
  cnv->cd(1);
  //
  TH1* dxodd = (TH1*)harr->At(np*10+kDTXodd);
  TH1* dxevn = (TH1*)harr->At(np*10+kDTXeven);
  TH1* dxoddS =(TH1*)harr->At(np*10+kDTXoddSPL);
  TH1* dxevnS =(TH1*)harr->At(np*10+kDTXevenSPL);
  double max = TMath::Max(dxodd->GetMaximum(),dxevn->GetMaximum());
  dxodd->SetMaximum(1.1*max);
  dxodd->GetXaxis()->SetTitle("#DeltaX, #mum");
  dxodd->SetTitle(Form("#DeltaX for clSize=%d",np));
  dxodd->Fit("gaus","","");
  dxevn->Fit("gaus","","sames");
  //
  dxoddS->Draw("sames");
  dxevnS->Draw("sames");
  //
  gPad->Modified();
  gPad->Update();
  SetStPadPos(dxodd,0.75,0.97,0.8,1., -1,dxodd->GetLineColor());
  SetStPadPos(dxevn,0.75,0.97,0.6,0.8, -1,dxevn->GetLineColor());
  SetStPadPos(dxoddS,0.75,0.97,0.4,0.6, -1,dxoddS->GetLineColor());
  SetStPadPos(dxevnS,0.75,0.97,0.2,0.4, -1,dxevnS->GetLineColor());
  //
  cnv->cd(2);
  TH1* dz  = (TH1*)harr->At(np*10+kDTZ);
  dz->SetTitle(Form("#DeltaZ for clSize=%d",np));
  dz->GetXaxis()->SetTitle("#DeltaZ, #mum");
  dz->Fit("gaus");
  TH1* dzS = (TH1*)harr->At(np*10+kDTZSPL);
  dz->Draw("sames");
  gPad->Modified();
  gPad->Update();
  SetStPadPos(dz,0.75,0.97,0.8,1., -1, dz->GetLineColor());
  SetStPadPos(dzS,0.75,0.97,0.5,0.7, -1, dzS->GetLineColor());
  gPad->Modified();
  gPad->Update();
  //
  cnv->cd();
  return cnv;
}


TPaveStats* GetStPad(TH1* hst)
{
  TList *lst = hst->GetListOfFunctions();
  if (!lst) return 0;
  int nf = lst->GetSize();
  for (int i=0;i<nf;i++) {
    TPaveStats *fnc = (TPaveStats*) lst->At(i);
    if (fnc->InheritsFrom("TPaveStats")) return fnc;
  }
  return 0;
  //
}


TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl, Int_t col)
{
  TPaveStats* pad = GetStPad(hst);
  if (!pad) return 0;
  pad->SetX1NDC( x1 );
  pad->SetX2NDC( x2 );
  pad->SetY1NDC( y1 );
  pad->SetY2NDC( y2 );
  if (stl>=0) pad->SetFillStyle(stl);
  if (col>=0) pad->SetTextColor(col);
  pad->SetFillColor(0);
  //
  gPad->Modified();
  return pad;
}
Commit	Line	Data
e7c0d8c7	1	TObjArray histoArr;
	2	enum {kNPixAll=0,kNPixSPL=1,kDR=0,kDTXodd,kDTXeven,kDTZ, kDTXoddSPL,kDTXevenSPL,kDTZSPL};
	3
	4	TPaveStats* GetStPad(TH1* hst);
	5	TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl=-1,Int_t col=-1);
	6	TCanvas* DrawNP(int np, TObjArray* harr=0, TCanvas* cnv=0);
	7	TH1* GetHistoClSize(int npix,int id,TObjArray* harr=0);
	8	void DrawReport(const char* psname, TObjArray* harr=0);
	9
4c1ef97c	10
	11	typedef struct {
	12	Int_t evID;
	13	Int_t volID;
	14	Int_t lrID;
	15	Int_t clID;
	16	Int_t nPix;
	17	Int_t nX;
	18	Int_t nZ;
	19	Float_t pt;
	20	Float_t eta;
	21	Float_t xyz[3];
	22	Float_t dX;
	23	Float_t dZ;
	24	Bool_t split;
	25	Bool_t prim;
	26	Int_t pdg;
	27	Int_t ntr;
	28	} clSumm;
	29
e7c0d8c7	30	void compClusHits(int nev=-1)
	31	{
	32	const int kSplit=0x1<<22;
	33	const int kSplCheck=0x1<<23;
	34	//
	35	gSystem->Load("libITSUpgradeBase");
	36	gSystem->Load("libITSUpgradeSim");
	37	gSystem->Load("libITSUpgradeRec");
	38	gROOT->SetStyle("Plain");
	39
	40	gAlice=NULL;
	41	AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
	42	runLoader->LoadgAlice();
	43
	44	gAlice = runLoader->GetAliRun();
	45
	46	runLoader->LoadHeader();
	47	runLoader->LoadKinematics();
	48	runLoader->LoadRecPoints();
	49	runLoader->LoadSDigits();
	50	runLoader->LoadHits();
	51
	52	AliLoader *dl = runLoader->GetDetectorLoader("ITS");
	53
	54	AliGeomManager::LoadGeometry("geometry.root");
	55	AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
	56	Int_t nLayers = gm->GetNLayers();
	57	AliITSUClusterPix::SetGeom(gm);
	58	//
	59	AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
	60	its->CreateClusterArrays();
	61	//
	62	Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
	63	//
	64	TTree * cluTree = 0x0;
	65	TTree *hitTree = 0x0;
	66	TClonesArray *hitList=new TClonesArray("AliITSUHit");
	67	//
	68	TObjArray arrMCTracks; // array of hit arrays for each particle
	69	//
	70	Float_t xyzClGloF[3];
	71	Double_t xyzClGlo[3],xyzClTr[3];
	72	Int_t labels[3];
	73	int nLab = 0;
	74	int nlr=its->GetNLayersActive();
	75	int ntotev = (Int_t)runLoader->GetNumberOfEvents();
	76	printf("N Events : %i \n",ntotev);
	77	if (nev>0) ntotev = TMath::Min(nev,ntotev);
4c1ef97c	78	//
	79	// output tree
	80	TFile* flOut = TFile::Open("clInfo.root","recreate");
	81	TTree* trOut = new TTree("clitsu","clitsu");
	82	clSumm cSum;
	83	trOut->Branch("evID", &cSum.evID ,"evID/I");
	84	trOut->Branch("volID",&cSum.volID,"volID/I");
	85	trOut->Branch("lrID", &cSum.lrID ,"lrID/I");
	86	trOut->Branch("clID", &cSum.clID ,"clID/I");
	87	trOut->Branch("nPix", &cSum.nPix ,"nPix/I");
	88	trOut->Branch("nX" , &cSum.nX ,"nX/I");
	89	trOut->Branch("nZ" , &cSum.nZ ,"nZ/I");
	90	trOut->Branch("pt" , &cSum.pt ,"pt/F");
	91	trOut->Branch("eta" ,&cSum.eta ,"eta/F");
	92	trOut->Branch("xyz", cSum.xyz, "xyz[3]/F");
	93	trOut->Branch("dX" , &cSum.dX ,"dX/F");
	94	trOut->Branch("dZ" , &cSum.dZ ,"dZ/F");
	95	trOut->Branch("split",&cSum.split,"split/O");
	96	trOut->Branch("prim", &cSum.prim, "prim/O");
	97	trOut->Branch("pdg", &cSum.pdg, "pdg/I");
	98	trOut->Branch("ntr", &cSum.ntr, "ntr/I");
e7c0d8c7	99	//
	100	for (Int_t iEvent = 0; iEvent < ntotev; iEvent++) {
	101	printf("\n Event %i \n",iEvent);
	102	runLoader->GetEvent(iEvent);
4c1ef97c	103	AliStack *stack = runLoader->Stack();
e7c0d8c7	104	cluTree=dl->TreeR();
	105	hitTree=dl->TreeH();
	106	hitTree->SetBranchAddress("ITS",&hitList);
	107	//
	108	// read clusters
	109	for (int ilr=nlr;ilr--;) {
	110	TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr));
	111	if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);}
	112	br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress());
	113	}
	114	cluTree->GetEntry(0);
	115	its->ProcessClusters();
	116	//
	117	// read hits
	118	for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
	119	hitTree->GetEntry(iEnt);
	120	int nh = hitList->GetEntries();
	121	for(Int_t iHit=0; iHit<hitList->GetEntries();iHit++){
	122	AliITSUHit pHit = (AliITSUHit)hitList->At(iHit);
	123	int mcID = pHit->GetTrack();
	124	TClonesArray* harr = arrMCTracks.GetEntriesFast()>mcID ? (TClonesArray*)arrMCTracks.At(mcID) : 0;
	125	if (!harr) {
	126	harr = new TClonesArray("AliITSUHit"); // 1st encounter of the MC track
	127	arrMCTracks.AddAtAndExpand(harr,mcID);
	128	}
	129	//
	130	new ( (harr)[harr->GetEntriesFast()] ) AliITSUHit(pHit);
	131	}
	132	}
	133	//
	134	// compare clusters and hits
	135	//
	136	printf(" tree entries: %d\n",cluTree->GetEntries());
	137	//
	138	for (int ilr=0;ilr<nlr;ilr++) {
	139	AliITSURecoLayer* lr = its->GetLayerActive(ilr);
	140	TClonesArray* clr = lr->GetClusters();
	141	int nClu = clr->GetEntries();
	142	printf("Layer %d : %d clusters\n",ilr,nClu);
	143	//
	144	for (int icl=0;icl<nClu;icl++) {
	145	AliITSUClusterPix cl = (AliITSUClusterPix)clr->At(icl);
	146	int modID = cl->GetVolumeId();
	147
	148	//------------ check if this is a split cluster
	149	int sInL = modID - gm->GetFirstModIndex(ilr);
	150	if (!cl->TestBit(kSplCheck)) {
	151	cl->SetBit(kSplCheck);
	152	// check if there is no other cluster with same label on this module
	153	AliITSURecoSens* sens = lr->GetSensor(sInL);
	154	int nclSn = sens->GetNClusters();
	155	int offs = sens->GetFirstClusterId();
	156	// printf("To check for %d (mod:%d) N=%d from %d\n",icl,modID,nclSn,offs);
	157	for (int ics=0;ics<nclSn;ics++) {
	158	AliITSUClusterPix* clusT = (AliITSUClusterPix*)lr->GetCluster(offs+ics); // access to clusters
	159	if (clusT==cl) continue;
	160	for (int ilb0=0;ilb0<3;ilb0++) {
	161	int lb0 = cl->GetLabel(ilb0); if (lb0<=-1) break;
	162	for (int ilb1=0;ilb1<3;ilb1++) {
	163	int lb1 = clusT->GetLabel(ilb1); if (lb1<=-1) break;
	164	if (lb1==lb0) {
	165	cl->SetBit(kSplit);
	166	clusT->SetBit(kSplit);
	167	/*
168	printf("Discard clusters of module %d:\n",modID);
169	cl->Print();
170	clusT->Print();
171	*/
172	break;
173	}
174	}
175	}
176	}
177	}
178	//------------
179	AliITSsegmentation* segm = gm->GetSegmentation(ilr);
180	//
181	cl->GetGlobalXYZ(xyzClGloF);
182	int clsize = cl->GetNPix();
183	for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
184	TGeoHMatrix* mat = gm->GetMatrixSens(modID);
185	if (!mat) {printf("failed to get matrix for module %d\n",cl->GetVolumeId());}
186	mat->MasterToLocal(xyzClGlo,xyzClTr);
187	//
188	int col,row;
189	segm->LocalToDet(xyzClTr[0],xyzClTr[2],row,col); // effective col/row
190	nLab = 0;
191	for (int il=0;il<3;il++) {
192	if (cl->GetLabel(il)>=0) labels[nLab++] = cl->GetLabel(il);
193	else break;
194	}
195	// find hit info
196	for (int il=0;il<nLab;il++) {
197	TClonesArray* htArr = (TClonesArray*)arrMCTracks.At(labels[il]);
198	if (!htArr) {printf("did not find MChits for label %d ",labels[il]); cl->Print(); continue;}
199	//
200	int nh = htArr->GetEntriesFast();
201	AliITSUHit *pHit=0;
202	for (int ih=nh;ih--;) {
203	AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih);
204	if (tHit->GetModule()!=modID) continue;
205	pHit = tHit;
206	break;
207	}
208	if (!pHit) {
209	printf("did not find MChit for label %d on module %d ",il,modID);
210	cl->Print();
211	htAtt->Print();
212	continue;
213	}
214	//
215	pHit->GetPositionG(xg1,yg1,zg1);
216	pHit->GetPositionG0(xg0,yg0,zg0,tg0);
217	//
218	double txyzH[3],gxyzH[3] = { (xg1+xg0)/2, (yg1+yg0)/2, (zg1+zg0)/2 };
219	mat->MasterToLocal(gxyzH,txyzH);
220	double rcl = TMath::Sqrt(xyzClTr[0]xyzClTr[0]+xyzClTr[1]xyzClTr[1]);
221	double rht = TMath::Sqrt(txyzH[0]txyzH[0]+txyzH[1]txyzH[1]);
222	//
223	GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
224	if (cl->TestBit(kSplit)) {
225	if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
226	else GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
227	GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
228	GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize);
229	}
230	if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
231	else GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
232	GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
233	GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize);
4c1ef97c	234	//
	235	cSum.evID = iEvent;
	236	cSum.volID = cl->GetVolumeId();
	237	cSum.lrID = ilr;
	238	cSum.clID = icl;
	239	cSum.nPix = cl->GetNPix();
	240	cSum.nX = cl->GetNx();
	241	cSum.nZ = cl->GetNz();
	242	cSum.split = cl->TestBit(kSplit);
	243	cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
	244	cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
	245	int label = cl->GetLabel(0);
	246	TParticle* part = 0;
	247	if (label>=0 && (part=stack->Particle(label)) ) {
	248	cSum.pdg = part->GetPdgCode();
	249	cSum.eta = part->Eta();
	250	cSum.pt = part->Pt();
	251	cSum.prim = stack->IsPhysicalPrimary(label);
	252	}
	253	cSum.ntr = 0;
	254	for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
	255	for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
	256	//
	257	trOut->Fill();
e7c0d8c7	258	/*
	259	if (clsize==5) {
	260	printf("\nL%d(%c) Mod%d, Cl:%d \| %+5.1f %+5.1f (%d/%d)\|H:%e %e %e \| C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
	261	modID,icl,(txyzH[0]-xyzClTr[0])1e4,(txyzH[2]-xyzClTr[2])1e4, row,col,
	262	gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
	263	cl->Print();
	264	pHit->Print();
	265	//
	266	double a0,b0,c0,a1,b1,c1,e0;
	267	pHit->GetPositionL0(a0,b0,c0,e0);
	268	pHit->GetPositionL(a1,b1,c1);
	269	float cloc[3];
	270	cl->GetLocalXYZ(cloc);
	271	printf("LocH: %e %e %e \| %e %e %e\n",a0,b0,c0,a1,b1,c1);
	272	printf("LocC: %e %e %e \| %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
	273	}
	274	*/
	275	//
	276	}
	277	}
	278	}
	279
	280	// layerClus.Clear();
	281	//
	282	arrMCTracks.Delete();
	283	}//event loop
	284	//
4c1ef97c	285	flOut->cd();
	286	trOut->Write();
	287	delete trOut;
	288	flOut->Close();
	289	flOut->Delete();
e7c0d8c7	290	DrawReport("clinfo.ps",&histoArr);
	291	//
	292	}
	293
	294	void DrawReport(const char* psname, TObjArray* harr)
	295	{
	296	gStyle->SetOptFit(1);
	297	if (!harr) harr = &histoArr;
	298	TCanvas* cnv = new TCanvas("cl","cl",900,600);
	299	//
	300	TString psnm1 = psname;
	301	if (psnm1.IsNull()) psnm1 = "clusters.ps";
	302	TString psnm0 = psnm1.Data();
	303	psnm0 += "[";
	304	TString psnm2 = psnm1.Data();
	305	psnm2 += "]";
	306	cnv->Print(psnm0.Data());
	307	//
	308	TH1* clall = GetHistoClSize(0,kNPixAll,harr);
	309	clall->SetLineColor(kRed);
	310	clall->Draw();
	311	TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr);
	312	clSpl->SetLineColor(kBlue);
	313	clSpl->Draw("sames");
	314	gPad->Modified();
	315	gPad->Update();
	316	SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor());
	317	SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor());
	318	gPad->Modified();
	319	gPad->Update();
	320	gPad->SetLogy(1);
	321	//
	322	cnv->cd();
	323	cnv->Print(psnm1.Data());
	324	//
	325	// plot cluster sized from 1 to 10
	326	for (int i=1;i<=10;i++) {
	327	if (clall->GetBinContent(clall->FindBin(i))<100) continue;
	328	DrawNP(i,harr,cnv);
	329	cnv->Print(psnm1.Data());
	330	}
	331	cnv->Print(psnm2.Data());
	332	}
	333
	334	//------------------------------------------
	335	TH1* GetHistoClSize(int npix,int id,TObjArray* harr)
	336	{
	337	// book histos
	338	TH1* h = 0;
	339	if (!harr) harr = &histoArr;
	340	//
	341	if (npix<1) {
	342	if (harr->GetEntriesFast()>=id && (h=(TH1*)harr->At(id))) return h;
	343	h = new TH1F("npixAll","npixAll",150,0.5,54.5);
	344	h->SetLineColor(kRed);
	345	harr->AddAtAndExpand(h, kNPixAll);
	346	//
	347	h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
	348	h->SetLineColor(kBlue);
	349	harr->AddAtAndExpand(h, kNPixSPL);
	350	//
	351	h = (TH1*)harr->At(id);
	352	if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
	353	return h;
354	}
355	//
356	int idh = npix*10+id;
357	if (harr->GetEntriesFast()>=idh && (h=(TH1*)harr->At(idh))) return h;
358	//
359	const int nbin=100;
360	const double kdiff=80;
361	// need to create set of histos
362	//
363	h = new TH1F(Form("dxy_npix%d",npix),Form("dr_npix%d",npix),nbin,-kdiff,kdiff);
364	harr->AddAtAndExpand(h, npix*10 + kDR);
365	//
366	h = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
367	h->SetLineColor(kRed);
368	harr->AddAtAndExpand(h, npix*10 + kDTXodd);
369	h = new TH1F(Form("dtxEVN_npix%d",npix),Form("dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
370	h->SetLineColor(kBlue);
371	harr->AddAtAndExpand(h, npix*10 + kDTXeven);
372	//
373	h = new TH1F(Form("dtz_npix%d",npix),Form("dtz_npix%d",npix),nbin,-kdiff,kdiff);
374	h->SetLineColor(kGreen);
375	harr->AddAtAndExpand(h, npix*10 + kDTZ);
376	//
377	//
378	h = new TH1F(Form("SPL_dtxODD_npix%d",npix),Form("SPL_dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
379	h->SetLineColor(kMagenta);
380	h->SetFillColor(kMagenta);
381	h->SetFillStyle(3001);
382	h->SetLineStyle(2);
383	harr->AddAtAndExpand(h, npix*10 + kDTXoddSPL);
384	h = new TH1F(Form("SPL_dtxEVN_npix%d",npix),Form("SPL_dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
385	h->SetLineColor(kCyan);
386	h->SetFillColor(kCyan);
387	h->SetFillStyle(3006);
388	h->SetLineStyle(2);
389	harr->AddAtAndExpand(h, npix*10 + kDTXevenSPL);
390	//
391	h = new TH1F(Form("SPL_dtz_npix%d",npix),Form("SPLdtz_npix%d",npix),nbin,-kdiff,kdiff);
392	harr->AddAtAndExpand(h, npix*10 + kDTZSPL);
393	//
394	h->SetLineColor(kGreen+2);
395	h->SetFillColor(kGreen+2);
396	h->SetLineStyle(2);
397	h->SetFillStyle(3001);
398	h = (TH1*)harr->At(idh);
399	if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
400	return h;
401	}
402
403
404
405	TCanvas* DrawNP(int np, TObjArray* harr, TCanvas* cnv)
406	{
407	if (!harr) harr = &histoArr;
408	if (!cnv) cnv = new TCanvas(Form("cnv%d",np),Form("cnv%d",np),900,700);
409	cnv->Clear();
410	cnv->Divide(2,1);
411	cnv->cd(1);
412	//
413	TH1* dxodd = (TH1)harr->At(np10+kDTXodd);
414	TH1* dxevn = (TH1)harr->At(np10+kDTXeven);
415	TH1* dxoddS =(TH1)harr->At(np10+kDTXoddSPL);
416	TH1* dxevnS =(TH1)harr->At(np10+kDTXevenSPL);
417	double max = TMath::Max(dxodd->GetMaximum(),dxevn->GetMaximum());
418	dxodd->SetMaximum(1.1*max);
419	dxodd->GetXaxis()->SetTitle("#DeltaX, #mum");
420	dxodd->SetTitle(Form("#DeltaX for clSize=%d",np));
421	dxodd->Fit("gaus","","");
422	dxevn->Fit("gaus","","sames");
423	//
424	dxoddS->Draw("sames");
425	dxevnS->Draw("sames");
426	//
427	gPad->Modified();
428	gPad->Update();
429	SetStPadPos(dxodd,0.75,0.97,0.8,1., -1,dxodd->GetLineColor());
430	SetStPadPos(dxevn,0.75,0.97,0.6,0.8, -1,dxevn->GetLineColor());
431	SetStPadPos(dxoddS,0.75,0.97,0.4,0.6, -1,dxoddS->GetLineColor());
432	SetStPadPos(dxevnS,0.75,0.97,0.2,0.4, -1,dxevnS->GetLineColor());
433	//
434	cnv->cd(2);
435	TH1* dz = (TH1)harr->At(np10+kDTZ);
436	dz->SetTitle(Form("#DeltaZ for clSize=%d",np));
437	dz->GetXaxis()->SetTitle("#DeltaZ, #mum");
438	dz->Fit("gaus");
439	TH1* dzS = (TH1)harr->At(np10+kDTZSPL);
440	dz->Draw("sames");
441	gPad->Modified();
442	gPad->Update();
443	SetStPadPos(dz,0.75,0.97,0.8,1., -1, dz->GetLineColor());
444	SetStPadPos(dzS,0.75,0.97,0.5,0.7, -1, dzS->GetLineColor());
445	gPad->Modified();
446	gPad->Update();
447	//
448	cnv->cd();
449	return cnv;
450	}
451
452
453	TPaveStats* GetStPad(TH1* hst)
454	{
455	TList *lst = hst->GetListOfFunctions();
456	if (!lst) return 0;
457	int nf = lst->GetSize();
458	for (int i=0;i<nf;i++) {
459	TPaveStats fnc = (TPaveStats) lst->At(i);
460	if (fnc->InheritsFrom("TPaveStats")) return fnc;
461	}
462	return 0;
463	//
464	}
465
466
467	TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl, Int_t col)
468	{
469	TPaveStats* pad = GetStPad(hst);
470	if (!pad) return 0;
471	pad->SetX1NDC( x1 );
472	pad->SetX2NDC( x2 );
473	pad->SetY1NDC( y1 );
474	pad->SetY2NDC( y2 );
475	if (stl>=0) pad->SetFillStyle(stl);
476	if (col>=0) pad->SetTextColor(col);
477	pad->SetFillColor(0);
478	//
479	gPad->Modified();
480	return pad;
481	}
482