--- /dev/null
+/*
+This macro loops over all the cluster (output of compClusHits.C) and creates a database containing the topologies of the cluster and infromation about:
+pattern in TBits format, coordinate of the centre of the pixel contining the COG wrt down-left corner of the osculating rectangle (fraction of the pixel pithc),
+distance between COG and centre of the pixel (fraction), number of fired pixels, number of colums, number of rows.
+*/
+
+#if !defined(__CInt_t__) || defined(__MAKECInt_t__)
+#include "AliESDEvent.h"
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TBits.h"
+#include "TArrayI.h"
+#include "TArrayF.h"
+#include "../ITS/UPGRADE/AliITSUClusterPix.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSUHit.h"
+#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
+#include "AliITSsegmentation.h"
+#include "AliGeomManager.h"
+
+#endif
+
+TTree* treeInp = 0;
+TTree* cluTree = 0;
+AliRunLoader* runLoader = 0;
+AliESDEvent* esd=0;
+
+AliITSUGeomTGeo* gm=0;
+Int_t nLayers = 0;
+AliITSURecoDet *its = 0;
+TObjArray patterns;
+TArrayF pattCount;
+Int_t nPatterns=0,nClusters=0;
+TString outDBFileName = "clusterTopology.root";
+
+void ProcessEvent(Int_t iev);
+void ProcChunk(const char* path);
+void BuildClTopDB(const char* inpData="AliESDs.root", const char* outDBFile="clusterTopology.root");
+void AccountTopology(TBits& patt);
+void CreateDB();
+
+void BuildClTopDB(const char* inpData, const char* outDBFile)
+{
+ //
+ // parse input list of single root file and process it
+ outDBFileName = outDBFile;
+ if (outDBFileName.IsNull()) outDBFileName = "clusterTopology.root";
+ //
+ AliGeomManager::LoadGeometry("geometry.root");
+ gm = new AliITSUGeomTGeo(kTRUE);
+ AliITSUClusterPix::SetGeom(gm);
+ nLayers = gm->GetNLayers();
+ its = new AliITSURecoDet(gm, "ITSInt_terface");
+ its->CreateClusterArrays();
+ esd = new AliESDEvent();
+ //
+ //
+ TString inpDtStr = inpData;
+ if (inpDtStr.EndsWith(".root") || inpDtStr.EndsWith(".zip#")) {
+ ProcChunk(inpDtStr.Data());
+ }
+ else {
+ ifstream inpf(inpData);
+ if (!inpf.good()) {
+ printf("Failed on input filename %s\n",inpData);
+ return;
+ }
+ //
+ inpDtStr.ReadLine(inpf);
+ while ( !inpDtStr.IsNull() ) {
+ inpDtStr = inpDtStr.Strip(TString::kBoth,' ');
+ if (inpDtStr.BeginsWith("//") || inpDtStr.BeginsWith("#")) {inpDtStr.ReadLine(inpf); continue;}
+ inpDtStr = inpDtStr.Strip(TString::kBoth,',');
+ inpDtStr = inpDtStr.Strip(TString::kBoth,'"');
+ ProcChunk(inpDtStr.Data());
+ inpDtStr.ReadLine(inpf);
+ }
+ }
+ //
+ CreateDB();
+ //
+}
+
+void CreateDB()
+{
+ // sort patterns according to frequencies and store them
+ printf("Storing %d patterns in %s\n",nPatterns, outDBFileName.Data());
+ if (nPatterns<1) {
+ printf("No patterns found\n");
+ return;
+ }
+ TArrayI sortIndex;
+ sortIndex.Set(nPatterns);
+ TMath::Sort(nPatterns,pattCount.GetArray(),sortIndex.GetArray());
+ //
+ //create a file.txt containing topologies and corresponding frequecies
+ ofstream top("topologies.txt");
+
+ TObjArray arrStore(nPatterns);
+ TVectorF arrFreq(nPatterns);
+ TVectorF arrzCOGPix(nPatterns); //It's the cordinate of the pixel containing along the rows direction
+ TVectorF arrxCOGPix(nPatterns); //It's the COG cordinate along the columns direction
+ TVectorF arrzCOGshift(nPatterns); //It's the distance between COG and the centre of the pixel containing COG
+ TVectorF arrxCOGshift(nPatterns);
+ TVectorF arrNpix(nPatterns);
+ TVectorF arrNcol(nPatterns);
+ TVectorF arrNrow(nPatterns);
+
+ TCanvas* cnv = new TCanvas("cnv","Freq",900,600);
+ TH1F* hfreq = new TH1F("hfreq","Pattern frequency",300/2,-0.5, 300/2-0.5);
+ hfreq->SetStats(0);
+ hfreq->GetXaxis()->SetTitle("pattern ID");
+ hfreq->SetFillColor(kRed);
+ hfreq->SetFillStyle(3005);
+ TLegend* leg = new TLegend(70, 300000, 140,450000 ,"","");
+ leg->AddEntry((TObject*)0, Form("Number of patterns: %d",nPatterns), "");
+ leg->AddEntry((TObject*)0, "Number of events: 50", "");
+ leg->SetFillColor(0);
+ leg->SetBorderSize(0);
+
+ for (Int_t i=0;i<nPatterns;i++) {
+ Int_t ind = sortIndex[i];
+ Float_t fr = arrFreq[i] = Float_t(pattCount[ind])/nClusters;
+ arrStore.AddAt(patterns[ind],i);
+ UInt_t idd = patterns[ind]->GetUniqueID();
+
+ hfreq->Fill(i,pattCount[ind]);
+
+ Int_t firedPixels=0;
+ Int_t tempxCOG=0;
+ Int_t tempzCOG=0;
+
+ Int_t rc = idd>>16;//It's the rows span
+ Int_t cl = idd&0xffff;//It's the columns span
+ top << Form("\nPatt %4d R:%d C: %d Freq:%f \n\n",i,rc,cl,fr);
+ for (Int_t ir=rc;ir--;) {
+ top << "|";
+ for (Int_t ic=0; ic<cl; ic++) {
+ top << Form("%c",((TBits*)patterns[ind])->TestBitNumber(ir*cl+ic) ? '+':' ');
+ if(((TBits*)patterns[ind])->TestBitNumber(ir*cl+ic)){
+ firedPixels++;
+ tempxCOG+=ir;
+ tempzCOG+=ic;
+ }
+ }
+ top << ("|\n");
+ }
+ Float_t xsh=Float_t((tempxCOG%firedPixels))/firedPixels; //distance between COG end centre of the pixel containing COG
+ Float_t zsh=Float_t((tempzCOG%firedPixels))/firedPixels;
+ tempxCOG/=firedPixels;
+ tempzCOG/=firedPixels;
+ if(xsh>0.5){
+ tempxCOG+=1;
+ xsh-=1;
+ }
+ if(zsh>0.5){
+ tempzCOG+=1;
+ zsh-=1;
+ }
+ Float_t xc=(Float_t) tempxCOG+0.5;
+ Float_t zc=(Float_t) tempzCOG+0.5;
+ arrxCOGPix[i]=xc;
+ arrxCOGshift[i]=xsh;
+ arrzCOGPix[i]=zc;
+ arrzCOGshift[i]=zsh;
+ arrNpix[i]=firedPixels;
+ arrNrow[i]=rc;
+ arrNcol[i]=cl;
+
+ top << "\nxCentre: " << xc << " + " << xsh
+ <<" zCentre: " << zc << " + " << zsh <<"\n\n...............................................\n";
+ }
+
+ cnv->cd();
+ hfreq->Draw();
+ leg->Draw();
+ cnv->Print("Pattern_freq.jpg");
+
+ top.close();
+
+ printf("\nWe have %d clusters corresponding to %d patterns!!!\nEnjoy!!\n\n", nClusters, nPatterns);
+
+ /*
+ //Controlling if frequencies have been stored correctly
+ ofstream a("frequency.txt");
+ a << setw(15) << "i" << setw(20) << "Frequency" << endl;
+ for(Int_t i=0; i<nPatterns; i++) a << setw(15) << i << setw(20) << arrFreq[i] << endl;
+ a.close();
+ */
+ //
+ TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
+ flDB->WriteObject(&arrStore,"TopDB","kSingleKey");
+ flDB->WriteObject(&arrFreq, "TopFreq","kSingleKey");
+ flDB->WriteObject(&arrxCOGPix,"xCOG","kSingleKey");
+ flDB->WriteObject(&arrxCOGshift,"xShift","kSingleKey");
+ flDB->WriteObject(&arrzCOGPix,"zCOG","kSingleKey");
+ flDB->WriteObject(&arrzCOGshift,"zShift","kSingleKey");
+ flDB->WriteObject(&arrNpix,"NPix","kSingleKey");
+ flDB->WriteObject(&arrNrow,"NRow","kSingleKey");
+ flDB->WriteObject(&arrNcol,"NCol","kSingleKey");
+ //
+ flDB->Close();
+ delete flDB;
+}
+
+void ProcChunk(const char* path)
+{
+ //
+ TString dir=path;
+ //
+ printf("Processing %s\n",dir.Data());
+ //
+ if (dir.EndsWith("AliESDs.root")) {
+ dir.ReplaceAll("AliESDs.root","");
+ }
+ //
+ //
+ runLoader = AliRunLoader::Open(Form("%sgalice.root",dir.Data()));
+ if (!runLoader) {
+ printf("galice not found\n");
+ return;
+ }
+ runLoader->LoadgAlice();
+ gAlice = runLoader->GetAliRun();
+ // runLoader->LoadHeader();
+ // runLoader->LoadKinematics();
+ runLoader->LoadRecPoints("ITS");
+ // ESD
+ TFile* esdFile = TFile::Open(Form("%sAliESDs.root",dir.Data()));
+ if (!esdFile || !esdFile->IsOpen()) {
+ printf("Error in opening ESD file\n");
+ // runLoader->UnloadKinematics();
+ // runLoader->UnloadHeader();
+ runLoader->UnloadgAlice();
+ delete runLoader;
+ return;
+ }
+
+ treeInp = (TTree*) esdFile->Get("esdTree");
+ if (!treeInp) {
+ printf("Error: no ESD tree found\n");
+ // runLoader->UnloadKinematics();
+ // runLoader->UnloadHeader();
+ runLoader->UnloadgAlice();
+ delete runLoader;
+ return;
+ }
+ esd->ReadFromTree(treeInp);
+ //
+ for (Int_t iEv=0;iEv<runLoader->GetNumberOfEvents(); iEv++) {
+ ProcessEvent(iEv);
+ }
+ runLoader->UnloadRecPoints("all");
+ // runLoader->UnloadKinematics();
+ // runLoader->UnloadHeader();
+ runLoader->UnloadgAlice();
+ delete runLoader;
+ runLoader = 0;
+ esdFile->Close();
+ delete esdFile;
+}
+
+
+void ProcessEvent(Int_t iEv)
+{
+ // process single event clusters
+ runLoader->GetEvent(iEv);
+ treeInp->GetEvent(iEv);
+ cluTree = runLoader->GetTreeR("ITS",kFALSE);
+ if (!cluTree) {
+ printf("Did not get cluster tree for event %d\n",iEv);
+ return;
+ }
+ //
+ for (Int_t ilr=0;ilr<nLayers;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();
+ //
+ TBits currTop;
+ //
+ for (Int_t ilr=0;ilr<nLayers;ilr++) {
+ AliITSURecoLayer* lr = its->GetLayerActive(ilr);
+ TClonesArray* clr = lr->GetClusters();
+ Int_t nClu = clr->GetEntries();
+ printf("Ev%d Lr:%d Ncl:%d\n",iEv,ilr,nClu);
+ // printf("Layer %d : %d clusters\n",ilr,nClu);
+ //
+ for (Int_t icl=0;icl<nClu;icl++) {
+ AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(icl);
+ currTop.Clear();
+ UShort_t rs = (UShort_t)cl->GetPatternRowSpan();
+ UShort_t cs = (UShort_t)cl->GetPatternColSpan();
+ for (Int_t ir=0;ir<rs;ir++)
+ for (Int_t ic=0;ic<cs;ic++)
+ if (cl->TestPixel(ir,ic)) currTop.SetBitNumber(ir*cs+ic);
+ //
+ currTop.SetUniqueID( (rs<<16) + (cs));
+ //
+ AccountTopology(currTop);
+ }
+ }
+ //
+}
+
+//___________________________________________
+void AccountTopology(TBits& patt)
+{
+ // account topology in the database
+ Bool_t newPatt = kTRUE;
+ nClusters++;
+ for (Int_t ip=0;ip<nPatterns;ip++) {
+ TBits* pattOld = (TBits*)patterns.At(ip);
+ if (*pattOld == patt && pattOld->GetUniqueID()==patt.GetUniqueID()) {
+ // require that the matrix size and the bit patterns are equal
+ newPatt = kFALSE;
+ pattCount[ip]++;
+ break;
+ }
+ }
+ //
+ if (newPatt) {
+ TBits* pt = new TBits(patt);
+ patterns.AddLast(pt);
+ if (pattCount.GetSize()<nPatterns+100) pattCount.Set(nPatterns+100);
+ pattCount[nPatterns++] = 1;
+ }
+ //
+}
\ No newline at end of file
--- /dev/null
+Macros from Luca Barioglio for claster patterns analysis
+Presented in
+https://indico.cern.ch/event/339513/contribution/3/material/slides/0.pdf
+- BuildClTopDB.C
+- compClusHits.C (modified with impact angles)
+- SecondcompClusHits.C (with DB Loading)
+- complete1.C (all the histograms)
+- errClus.C (just 1D histograms)
+- groupdef.C (grouping)
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TBits.h"
+#include "../ITS/UPGRADE/AliITSUClusterPix.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSUHit.h"
+#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
+#include "AliITSsegmentation.h"
+#include "AliGeomManager.h"
+#include "AliStack.h"
+#include "AliLoader.h"
+#include "AliCDBManager.h"
+
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TGeoMatrix.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TPaveStats.h"
+#include "TClonesArray.h"
+
+#endif
+
+TObjArray *pattDB=0;
+TVectorF* pattFR=0;
+TVectorF* xCentrePix=0;
+TVectorF* zCentrePix=0;
+TVectorF* xCentreShift=0;
+TVectorF* zCentreShift=0;
+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);
+void LoadDB(const char* fname);
+
+
+typedef struct {
+ Int_t evID;
+ Int_t volID;
+ Int_t lrID;
+ Int_t clID;
+ Int_t nPix;
+ Int_t nX;
+ Int_t nZ;
+ Int_t q;
+ Float_t pt;
+ Float_t eta;
+ Float_t phi;
+ Float_t xyz[3];
+ Float_t dX;
+ Float_t dY;
+ Float_t dZ;
+ Bool_t split;
+ Bool_t prim;
+ Int_t pdg;
+ Int_t ntr;
+ Float_t alpha; // alpha is the angle in y-radius plane in local frame
+ Float_t beta; // beta is the angle in xz plane, taken from z axis, growing counterclockwise
+ Int_t nRowPatt;
+ Int_t nColPatt;
+ Int_t pattID;
+ Float_t freq;
+ Float_t xCen;
+ Float_t zCen;
+ Float_t zShift;
+ Float_t xShift;
+} clSumm;
+
+TBits clTop;
+
+void SecondcompClusHits(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");
+
+ AliCDBManager* man = AliCDBManager::Instance();
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetSpecificStorage("GRP/GRP/Data",
+ Form("local://%s",gSystem->pwd()));
+ man->SetSpecificStorage("ITS/Align/Data",
+ Form("local://%s",gSystem->pwd()));
+ man->SetSpecificStorage("ITS/Calib/RecoParam",
+ Form("local://%s",gSystem->pwd()));
+ man->SetRun(0);
+
+
+ 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");
+ TObjArray algITS;
+ AliGeomManager::LoadAlignObjsFromCDBSingleDet("ITS",algITS);
+ AliGeomManager::ApplyAlignObjsToGeom(algITS);
+ //
+ AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+ AliITSUClusterPix::SetGeom(gm);
+ //
+ AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
+ its->CreateClusterArrays();
+ //
+ Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
+ Double_t xExit,yExit,zExit,xEnt,yEnt,zEnt,tof1;
+ //
+ 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);
+ //
+ // load topology database
+ LoadDB("clusterTopology.root");
+
+ Int_t nPatterns = pattDB->GetEntriesFast();
+
+ // 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("q" , &cSum.q ,"q/I");
+ trOut->Branch("pt" , &cSum.pt ,"pt/F");
+ trOut->Branch("eta" ,&cSum.eta ,"eta/F");
+ trOut->Branch("phi" , &cSum.phi ,"phi/F");
+ trOut->Branch("xyz", cSum.xyz, "xyz[3]/F");
+ trOut->Branch("dX" , &cSum.dX ,"dX/F");
+ trOut->Branch("dY" , &cSum.dY ,"dY/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");
+ trOut->Branch("alpha", &cSum.alpha, "alpha/F");
+ trOut->Branch("beta", &cSum.beta, "beta/F");
+ trOut->Branch("nRowPatt", &cSum.nRowPatt, "nRowPatt/I");
+ trOut->Branch("nColPatt", &cSum.nColPatt, "nColPatt/I");
+ trOut->Branch("pattID", &cSum.pattID,"pattID/I");
+ trOut->Branch("freq", &cSum.freq, "freq/F");
+ trOut->Branch("xCen", &cSum.xCen, "xCen/F");
+ trOut->Branch("zCen", &cSum.zCen, "zCen/F");
+ trOut->Branch("xShift", &cSum.xShift, "xShift/F");
+ trOut->Branch("zShift", &cSum.zShift, "zShift/F");
+ //
+ 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<nh;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: %lld\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();
+
+ //Compare clusters with DB
+
+ clTop.Clear();
+ Int_t tempPattID=-1;
+ Float_t tempFreq=0;
+ Float_t tempxCen=0;
+ Float_t tempzCen=0;
+ Float_t tempxShift=0;
+ Float_t tempzShift=0;
+ UShort_t rs = cl->GetPatternRowSpan();
+ UShort_t cs = cl->GetPatternColSpan();
+
+ for(int ir=0; ir<rs; ir++)
+ for(int ic=0; ic<cs; ic++)
+ if(cl->TestPixel(ir,ic)) clTop.SetBitNumber(ir*cs+ic);//Set cluster pattern
+ clTop.SetUniqueID((rs<<16)+cs);//Set cluster ID
+
+ for(int iDB=0; iDB<nPatterns; iDB++){
+ TBits* refPatt = (TBits*) pattDB->At(iDB);
+ if (*refPatt == clTop && refPatt->GetUniqueID()== clTop.GetUniqueID()){
+ tempPattID=iDB;
+ tempFreq=(*pattFR)[iDB];
+ tempxCen=(*xCentrePix)[iDB];
+ tempzCen=(*zCentrePix)[iDB];
+ tempxShift=(*xCentreShift)[iDB];
+ tempzShift=(*zCentreShift)[iDB];
+ break;
+ }
+ }
+
+ //------------ check if this is a split cluster
+ if (!cl->TestBit(kSplCheck)) {
+ cl->SetBit(kSplCheck);
+ // check if there is no other cluster with same label on this module
+ AliITSURecoSens* sens = lr->GetSensorFromID(modID);
+ 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;
+ }
+ }
+ }
+ }
+ }
+ //------------
+ const AliITSsegmentation* segm = gm->GetSegmentation(ilr);
+ //
+ cl->GetGlobalXYZ(xyzClGloF);
+ int clsize = cl->GetNPix();
+ for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
+ const 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;
+ double dst2Max = 1e33;
+ for (int ih=nh;ih--;) {
+ AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih);
+ if (tHit->GetChip()!=modID) continue;
+ tHit->GetPositionG(xg1,yg1,zg1);
+ tHit->GetPositionG0(xg0,yg0,zg0,tg0);
+ double gxyzHDif[3] = { (xg1+xg0)/2 - xyzClGlo[0], (yg1+yg0)/2 - xyzClGlo[1], (zg1+zg0)/2 - xyzClGlo[2] };
+ double dst2 = gxyzHDif[0]*gxyzHDif[0] + gxyzHDif[1]*gxyzHDif[1] + gxyzHDif[2]*gxyzHDif[2];
+ if (dst2<dst2Max) {
+ pHit = tHit;
+ dst2Max = dst2;
+ }
+ }
+ if (!pHit) {
+ printf("did not find MChit for label %d on module %d ",il,modID);
+ cl->Print();
+ htArr->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]);
+ //
+ //
+ //Angles determination
+
+ pHit->GetPositionL(xExit,yExit,zExit);
+ pHit->GetPositionL0(xEnt,yEnt,zEnt,tof1);
+
+ Double_t dirHit[3]={(xExit-xEnt),(yExit-yEnt),(zExit-zEnt)};
+
+ Double_t alpha1 = TMath::ACos(TMath::Abs(dirHit[1])/TMath::Sqrt(dirHit[0]*dirHit[0]+dirHit[1]*dirHit[1]+dirHit[2]*dirHit[2])); //Polar Angle
+ Float_t alpha2 = (Float_t) alpha1; //convert to float
+ cSum.alpha = alpha2;
+
+ Double_t beta1;
+ beta1 = TMath::ATan2(TMath::Abs(dirHit[0]),TMath::Abs(dirHit[2])); //Azimuthal angle, values from -Pi to Pi
+ Float_t beta2 = (Float_t) beta1;
+ cSum.beta = beta2;
+
+ 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.q = cl->GetQ();
+ cSum.split = cl->TestBit(kSplit);
+ cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
+ cSum.dY = (txyzH[1]-xyzClTr[1])*1e4;
+ cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
+ cSum.nRowPatt = cl-> GetPatternRowSpan();
+ cSum.nColPatt = cl-> GetPatternColSpan();
+ cSum.pattID = tempPattID;
+ cSum.freq = tempFreq;
+ cSum.xCen = tempxCen;
+ cSum.zCen = tempzCen;
+ cSum.xShift = tempxShift;
+ cSum.zShift = tempzShift;
+
+ 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.phi = part->Phi();
+ 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->SetDirectory(0);
+ h->SetLineColor(kRed);
+ harr->AddAtAndExpand(h, kNPixAll);
+ //
+ h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
+ h->SetLineColor(kBlue);
+ h->SetDirectory(0);
+ harr->AddAtAndExpand(h, kNPixSPL);
+ //
+ h = (TH1*)harr->At(id);
+ if (!h) {printf("Unknown histo id=%d\n",id); 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);
+ h->SetDirectory(0);
+ harr->AddAtAndExpand(h, npix*10 + kDR);
+ //
+ h = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
+ h->SetDirectory(0);
+ 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->SetDirectory(0);
+ 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);
+ h->SetDirectory(0);
+ 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);
+ h->SetDirectory(0);
+
+ 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);
+ h->SetDirectory(0);
+ 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->SetDirectory(0);
+ //
+ 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;
+}
+
+void LoadDB(const char* fname)
+{
+ // load database
+ TFile* fl = TFile::Open(fname);
+ pattDB = (TObjArray*)fl->Get("TopDB");
+ pattFR = (TVectorF*)fl->Get("TopFreq");
+ xCentrePix =(TVectorF*)fl->Get("xCOG");
+ zCentrePix =(TVectorF*)fl->Get("zCOG");
+ xCentreShift =(TVectorF*)fl->Get("xShift");
+ zCentreShift =(TVectorF*)fl->Get("zShift");
+
+}
+
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "../ITS/UPGRADE/AliITSUClusterPix.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSUHit.h"
+#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
+#include "AliITSsegmentation.h"
+#include "AliGeomManager.h"
+#include "AliStack.h"
+#include "AliLoader.h"
+#include "AliCDBManager.h"
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TGeoMatrix.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TPaveStats.h"
+#include "TClonesArray.h"
+
+#endif
+
+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;
+ Int_t q;
+ Float_t pt;
+ Float_t eta;
+ Float_t phi;
+ Float_t xyz[3];
+ Float_t dX;
+ Float_t dY;
+ Float_t dZ;
+ Bool_t split;
+ Bool_t prim;
+ Int_t pdg;
+ Int_t ntr;
+ Float_t alpha; // alpha is the angle in y-radius plane in local frame
+ Float_t beta; // beta is the angle in xz plane, taken from z axis, growing counterclockwise
+ Int_t nRowPatt;
+ Int_t nColPatt;
+} 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");
+
+ AliCDBManager* man = AliCDBManager::Instance();
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetSpecificStorage("GRP/GRP/Data",
+ Form("local://%s",gSystem->pwd()));
+ man->SetSpecificStorage("ITS/Align/Data",
+ Form("local://%s",gSystem->pwd()));
+ man->SetSpecificStorage("ITS/Calib/RecoParam",
+ Form("local://%s",gSystem->pwd()));
+ man->SetRun(0);
+
+ 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");
+ TObjArray algITS;
+ AliGeomManager::LoadAlignObjsFromCDBSingleDet("ITS",algITS);
+ AliGeomManager::ApplyAlignObjsToGeom(algITS);
+ //
+ AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+ AliITSUClusterPix::SetGeom(gm);
+ //
+ AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
+ its->CreateClusterArrays();
+ //
+ Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
+ Double_t xExit,yExit,zExit,xEnt,yEnt,zEnt,tof1;
+
+ //
+ 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("q" , &cSum.q ,"q/I");
+ trOut->Branch("pt" , &cSum.pt ,"pt/F");
+ trOut->Branch("eta" ,&cSum.eta ,"eta/F");
+ trOut->Branch("phi" , &cSum.phi ,"phi/F");
+ trOut->Branch("xyz", cSum.xyz, "xyz[3]/F");
+ trOut->Branch("dX" , &cSum.dX ,"dX/F");
+ trOut->Branch("dY" , &cSum.dY ,"dY/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");
+ trOut->Branch("alpha", &cSum.alpha, "alpha/F");
+ trOut->Branch("beta", &cSum.beta, "beta/F");
+ trOut->Branch("nRowPatt", &cSum.nRowPatt, "nRowPatt/I");
+ trOut->Branch("nColPatt", &cSum.nColPatt, "nColPatt/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<nh;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: %lld\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->GetFirstChipIndex(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;
+ }
+ }
+ }
+ }
+ }
+ //------------
+ const AliITSsegmentation* segm = gm->GetSegmentation(ilr);
+ //
+ cl->GetGlobalXYZ(xyzClGloF);
+ int clsize = cl->GetNPix();
+ for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
+ const 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->GetChip()!=modID) continue;
+ pHit = tHit;
+ break;
+ }
+ if (!pHit) {
+ printf("did not find MChit for label %d on module %d ",il,modID);
+ cl->Print();
+ htArr->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]);
+ //
+ //Angles determination
+
+ pHit->GetPositionL(xExit,yExit,zExit);
+ pHit->GetPositionL0(xEnt,yEnt,zEnt,tof1);
+
+ Double_t dirHit[3]={(xExit-xEnt),(yExit-yEnt),(zExit-zEnt)};
+
+ /*double PG[3] = {(double)pHit->GetPXG(), (double)pHit->GetPYG(), (double)pHit->GetPZG()}; //Momentum at hit-point in Global Frame
+ double PL[3];
+ if (TMath::Abs(PG[0])<10e-7 && TMath::Abs(PG[1])<10e-7) {
+ pHit->Dump();
+ int lb = pHit->GetTrack();
+ stack->Particle(lb)->Print();
+ continue;
+ }
+ mat->MasterToLocalVect(PG,PL); //Momentum in local Frame
+ //printf(">> %e %e %e %e %e %e\n",PG[0],PL[0],PG[1],PL[1],PG[2],PL[2]);*/
+
+ Double_t alpha1 = TMath::ACos(TMath::Abs(dirHit[1])/TMath::Sqrt(dirHit[0]*dirHit[0]+dirHit[1]*dirHit[1]+dirHit[2]*dirHit[2])); //Polar Angle
+ Float_t alpha2 = (Float_t) alpha1; //convert to float
+ cSum.alpha = alpha2;
+
+ Double_t beta1;
+ beta1 = TMath::ATan2(dirHit[0],dirHit[2]); //Azimuthal angle, values from -Pi to Pi
+ Float_t beta2 = (Float_t) beta1;
+ cSum.beta = beta2;
+
+ 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.q = cl->GetQ();
+ cSum.split = cl->TestBit(kSplit);
+ cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
+ cSum.dY = (txyzH[1]-xyzClTr[1])*1e4;
+ cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
+ cSum.nRowPatt = cl-> GetPatternRowSpan();
+ cSum.nColPatt = cl-> GetPatternColSpan();
+
+ 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.phi = part->Phi();
+ 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->SetDirectory(0);
+ h->SetLineColor(kRed);
+ harr->AddAtAndExpand(h, kNPixAll);
+ //
+ h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
+ h->SetLineColor(kBlue);
+ h->SetDirectory(0);
+ harr->AddAtAndExpand(h, kNPixSPL);
+ //
+ h = (TH1*)harr->At(id);
+ if (!h) {printf("Unknown histo id=%d\n",id); 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);
+ h->SetDirectory(0);
+ harr->AddAtAndExpand(h, npix*10 + kDR);
+ //
+ h = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
+ h->SetDirectory(0);
+ 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->SetDirectory(0);
+ 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);
+ h->SetDirectory(0);
+ 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);
+ h->SetDirectory(0);
+
+ 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);
+ h->SetDirectory(0);
+ 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->SetDirectory(0);
+ //
+ 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;
+}
\ No newline at end of file
--- /dev/null
+/*
+This macro creates for each pattern four 2D histograms: distance between MChit and COG vs impact angle. This information are taken
+from the output of compClusHits.C. Information about the pattern are taken from the database generated by BuildClTop.C.
+The macro fits these histograms with gaussians along y axis, in order to study the dependence of dZ and dX wrt
+the imact angles. Two 1D histograms with dZ and dX are generated as a projection on to the y axis and are used to calculate mean value-sigma
+for each pattern.
+The output is a PDF with 4 2D histograms (mean-value of dZ(dX) vs alpha, sigma of dZ(dX) vs alpha,mean-value of dZ(dX) vs beta,
+sigma of dZ(dX) vs beta) and 2 histograms, dZ and dX, for each pattern, a drawing of the pattern and the information from the fit and the DB and
+information from DB and fits.
+You can set the range of patterns you want to study by setting LowerLimit/UpperLimit. If the range coincide with the the whole ragne of patterns,
+the DB is updated with the information from these fits.
+It is possible to loop over few clusters instead of all the cluster setting ntotCl: usefull to check the macro. (-1 means all the clusters).
+It is possilble toprint a PDF with the 2D histograms, pattern drawing and info: setting number2d histos you decide the number of pattern for which
+you want 2D histograms.
+*/
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TF1.h"
+#include "TArrayI.h"
+#include "TArrayF.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TPad.h"
+#include "TLatex.h"
+#include "TBits.h"
+#include "TStopwatch.h"
+#include "TLine.h"
+#include "TMarker.h"
+#include "../ITS/UPGRADE/AliITSUClusterPix.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSUHit.h"
+#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
+#include "AliITSsegmentation.h"
+#include "AliGeomManager.h"
+#include "TROOT.h"
+#include "TStyle.h"
+
+#endif
+
+TObjArray* pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
+TVectorF* pattFR=0; //frequency of the pattern in the database
+TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
+TVectorF* zCentrePix=0;
+TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
+TVectorF* zCentreShift=0;
+TVectorF* NPix=0; //Number of fired pixels
+TVectorF* NRow=0; //Number of rows of the pattern
+TVectorF* NCol=0; //Number of columns of the pattern
+
+TArrayF DeltaZmean; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TArrayF DeltaZmeanErr;
+TArrayF DeltaXmean; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TArrayF DeltaXmeanErr;
+TArrayF DeltaZsigma; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TArrayF DeltaZsigmaErr;
+TArrayF DeltaXsigma; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TArrayF DeltaXsigmaErr;
+TArrayF Chi2x; //Chi2 of the gaussian fit over x
+TArrayF Chi2z; //Chi2 of the aussian fut over z
+TArrayI NDFx;
+TArrayI NDFz;
+
+TObjArray histoArr;
+
+TStopwatch timer;
+TStopwatch timer1;
+
+Float_t pitchX; //pitch of the pixel in x direction
+Float_t pitchZ; //pitch of the pixel in z direction
+
+enum{kDXAlpha=0, kDZAlpha=1, kDXBeta=2, kDZBeta=3};//to select the kind of TH2 we want
+
+void LoadDB(const char* fname);
+void CreateHistos(TObjArray* harr);
+TH2* GetHisto(Int_t pattID, Int_t kind, TObjArray* harr);
+void FitHistos(Int_t pattID, Float_t Freq, Int_t nPix, Int_t nRow, Int_t nCol,
+ TObjArray* harr, TCanvas* canvas, TCanvas* canvasbis);
+void UpDateDB(const char* outDBFile);
+
+typedef struct {
+ Int_t evID;
+ Int_t volID;
+ Int_t lrID;
+ Int_t clID;
+ Int_t nPix;
+ Int_t nX;
+ Int_t nZ;
+ Int_t q;
+ Float_t pt;
+ Float_t eta;
+ Float_t phi;
+ Float_t xyz[3];
+ Float_t dX;
+ Float_t dY;
+ Float_t dZ;
+ Bool_t split;
+ Bool_t prim;
+ Int_t pdg;
+ Int_t ntr;
+ Float_t alpha; // alpha is the angle in y-radius plane in local frame
+ Float_t beta; // beta is the angle in xz plane, taken from z axis, growing counterclockwise
+ Int_t nRowPatt;
+ Int_t nColPatt;
+ Int_t pattID;
+ Float_t freq;
+ Float_t xCen;
+ Float_t zCen;
+ Float_t zShift;
+ Float_t xShift;
+} clSumm;
+
+static clSumm Summ;
+
+Int_t nPatterns=0;
+Int_t LowerLimit=0;//set the first cluster of the range of interest, INCLUDED (to update DB LowerLimit=0)
+Int_t UpperLimit=4000; // EXCLUDED; -1 means number of Patterns. In this case the database is updated with the information about DeltaxMean, DeltaXerr, etc.
+Int_t nStudied=0;
+Int_t ntotCl=-1;//Set -1 to loop over all the clusters, otherwise set the number of clusters
+Int_t number2dhisto=10;//Set the number pattern from 0 to number2dhisto-1 for which you want to print the TH2 histos (dZ or dX vs alpha or beta)
+
+void complete1(){
+
+ //importing data
+
+ LoadDB("clusterTopology.root");
+
+ AliGeomManager::LoadGeometry("geometry.root");
+ AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+ AliITSUClusterPix::SetGeom(gm);
+ const AliITSsegmentation* segm = gm->GetSegmentation(0);
+ pitchX = segm->Dpx(0)*10000; //pitch in X in micron
+ printf("pitchX: %f\n",pitchX);
+ pitchZ = segm->Dpz(0)*10000; //pitch in Z in micron
+ printf("pitchX: %f\n",pitchX);
+ delete gm;
+
+ nPatterns = (pattDB->GetEntriesFast());
+
+ if(UpperLimit==-1){
+ UpperLimit=nPatterns;
+ }
+
+ nStudied=UpperLimit-LowerLimit+1;
+
+ printf("\n\n nPatterns %d\n\n",nPatterns);
+
+ DeltaZmean.Set(nPatterns+100);
+ DeltaZmeanErr.Set(nPatterns+100);
+ DeltaXmean.Set(nPatterns+100);
+ DeltaXmeanErr.Set(nPatterns+100);
+ DeltaZsigma.Set(nPatterns+100);
+ DeltaZsigmaErr.Set(nPatterns+100);
+ DeltaXsigma.Set(nPatterns+100);
+ DeltaXsigmaErr.Set(nPatterns+100);
+ Chi2x.Set(nPatterns+100);
+ Chi2z.Set(nPatterns+100);
+ NDFx.Set(nPatterns+100);
+ NDFz.Set(nPatterns+100);
+
+ printf("Loading input tree... ");
+
+ TFile *input = new TFile("clInfo.root","read");
+ TTree *clitsu = (TTree*) input->Get("clitsu");
+
+ if(ntotCl==-1){
+ ntotCl = (Int_t) clitsu->GetEntries();
+ }
+ TBranch *brdZ = clitsu->GetBranch("dZ");
+ brdZ->SetAddress(&Summ.dZ);
+ TBranch *brdX = clitsu->GetBranch("dX");
+ brdX->SetAddress(&Summ.dX);
+ TBranch *brpattID = clitsu->GetBranch("pattID");
+ brpattID->SetAddress(&Summ.pattID);
+ TBranch *brfreq = clitsu->GetBranch("freq");
+ brfreq->SetAddress(&Summ.freq);
+ TBranch *bralpha = clitsu->GetBranch("alpha");
+ bralpha->SetAddress(&Summ.alpha);
+ TBranch *brbeta = clitsu->GetBranch("beta");
+ brbeta->SetAddress(&Summ.beta);
+ printf("done!!\n\n");
+
+ CreateHistos(&histoArr);
+ printf("Starting to process clusters\n\n");
+ timer.Start();
+ for(Int_t iCl=0; iCl<ntotCl; iCl++){
+ clitsu->GetEntry(iCl);
+ Int_t ID = Summ.pattID;
+ printf("Processing cluster %d... ",iCl);
+ if(ID<LowerLimit || ID>=UpperLimit) {
+ printf("skipped!\n\n");
+ continue;
+ }
+ GetHisto(ID,kDXAlpha,&histoArr)->Fill(Summ.alpha,Summ.dX);
+ GetHisto(ID,kDZAlpha,&histoArr)->Fill(Summ.alpha,Summ.dZ);
+ GetHisto(ID,kDXBeta,&histoArr)->Fill(Summ.beta,Summ.dX);
+ GetHisto(ID,kDZBeta,&histoArr)->Fill(Summ.beta,Summ.dZ);
+ printf("done!\n\n");
+ }
+ printf("All clusters processed!\n");
+ timer.Stop();
+ timer.Print();
+
+ TCanvas* cnv1 = new TCanvas("cnv1","cnv1",900,700);
+ TCanvas* cnv2 = new TCanvas("cnv2","cnv2",900,700);
+
+ printf("Printing PDF...\n\n");
+ timer1.Start();
+
+ cnv1->Print("angles.pdf[");
+ cnv2->Print("angles2D.pdf[");
+
+ for(Int_t i=LowerLimit; i<UpperLimit; i++){
+
+ FitHistos(i,(*pattFR)[i],(*NPix)[i],(*NRow)[i],(*NCol)[i],&histoArr, cnv1, cnv2);
+ }
+
+ cnv1->Print("angles.pdf]");
+ cnv2->Print("angles2D.pdf]");
+ timer1.Stop();
+ printf("\n\ndone!!\n\n");
+ timer1.Print();
+
+ delete cnv1;
+ delete cnv2;
+
+ if(LowerLimit==0 && UpperLimit==nPatterns){
+ UpDateDB("clusterTopology2.root");
+ }
+
+}
+
+void LoadDB(const char* fname){ //load the data base built by BuildClTopDB.C
+
+ // load database
+ TFile* fl = TFile::Open(fname);
+ if(!fl){printf("Could not find %s",fname); exit(1);}
+ pattDB = (TObjArray*)fl->Get("TopDB");
+ pattFR = (TVectorF*)fl->Get("TopFreq");
+ xCentrePix =(TVectorF*)fl->Get("xCOG");
+ zCentrePix =(TVectorF*)fl->Get("zCOG");
+ xCentreShift =(TVectorF*)fl->Get("xShift");
+ zCentreShift =(TVectorF*)fl->Get("zShift");
+ NPix =(TVectorF*)fl->Get("NPix");
+ NCol =(TVectorF*)fl->Get("NCol");
+ NRow =(TVectorF*)fl->Get("NRow");
+}
+
+void CreateHistos(TObjArray* harr){ //creates four 2D histograms (dX or dZ vs alpha or beta) for each pattern included in the study range
+
+ printf("Creating histograms... ");
+
+ //for each pattern we crate 4 2D histograms corresponding the cathegories within the enumeration
+ if(!harr) harr=&histoArr;
+
+ for(Int_t i=0; i<nStudied; i++){
+
+ TH2* h0 = new TH2F(Form("hXalpha%d", i),
+ Form("#DeltaX vs #alpha for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
+ h0->SetDirectory(0);
+ harr->AddAtAndExpand(h0, i*4);//This histogram is the first of the 4-histos series
+
+ TH2* h1 = new TH2F(Form("hZalpha%d", i),
+ Form("#DeltaZ vs #alpha for pattID %d", i),10, -0.1, 1.1* TMath::Pi()/2,100,-30,30);
+ h1->SetDirectory(0);
+ harr->AddAtAndExpand(h1, i*4+1);//This histogram is the second of the 4-histos series and so on...
+
+ TH2* h2 = new TH2F(Form("hXbeta%d", i),
+ Form("#DeltaX vs #beta for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
+ h2->SetDirectory(0);
+ harr->AddAtAndExpand(h2, i*4+2);
+
+ TH2* h3 = new TH2F(Form("hZbeta%d", i),
+ Form("#DeltaZ vs #beta for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
+ h3->SetDirectory(0);
+ harr->AddAtAndExpand(h3, i*4+3);
+
+ }
+
+ printf("done!!\n\n");
+}
+
+TH2* GetHisto(Int_t pattID, Int_t kind, TObjArray* harr){
+
+ if(!harr) harr=&histoArr ;
+ TH2* h=0;
+ h=(TH2*)harr->At((pattID-LowerLimit)*4+kind);
+ if (!h) {
+ printf("Unknown histo id=%d\n",pattID);
+ exit(1);
+ }
+ return h;
+
+}
+
+void FitHistos(Int_t pattID, Float_t Freq, Int_t nPix, Int_t nRow,
+ Int_t nCol, TObjArray* harr, TCanvas* canvas, TCanvas* canvasbis){ //Fit the 2D histograms with gaussian along the Y and creates two 1D histograms,
+ // with dX and dZ, bay a projection to the y axis and the fit with gaussian
+
+ static TF1* gs = new TF1("gs","gaus",-50,50);
+
+ if(!harr) harr=&histoArr;
+
+ //Creating the pads for the output
+
+ canvas->Clear();
+ canvas->cd();
+
+ Float_t width = (1 - canvas->GetRightMargin() - canvas->GetLeftMargin() )/3;
+
+ TPad* titlepad = new TPad("titlepad","",canvas->GetLeftMargin(),0.825,1-canvas->GetRightMargin(),0.95);
+ titlepad->SetFillColor(kYellow);
+ titlepad->Draw();
+ TPad* pad1 = new TPad("pad1","",canvas->GetLeftMargin(),(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),canvas->GetLeftMargin()+width,0.8);
+ pad1->Draw();
+ TPad* pad2 = new TPad("pad2","",canvas->GetLeftMargin()+width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),canvas->GetLeftMargin()+2*width,0.8);
+ pad2->Draw();
+ TPad* pad3 = new TPad("pad3","",canvas->GetLeftMargin()+2*width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),1-canvas->GetRightMargin(),0.8);
+ pad3->Draw();
+ TPad* pad4 = new TPad("pad3","",canvas->GetLeftMargin(),canvas->GetBottomMargin(),canvas->GetLeftMargin()+width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
+ pad4->Draw();
+ TPad* pad5 = new TPad("pad5","",canvas->GetLeftMargin()+width,canvas->GetBottomMargin(),canvas->GetLeftMargin()+2*width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
+ pad5->Draw();
+ TPad* pad6 = new TPad("pad6","",canvas->GetLeftMargin()+2*width,canvas->GetBottomMargin(),1-canvas->GetRightMargin(),(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
+ pad6->Draw();
+ TPad* padPattern = new TPad("padPattern","",0.75,0.825,0.9,0.95);
+ padPattern->Draw();
+
+ //...................................................................................................
+
+ // dX vs alpha
+
+ TObjArray harrFXA;
+ TH2* hXA = (TH2*)harr->At((pattID-LowerLimit)*4);
+ hXA->GetYaxis()->SetTitle("#DeltaX, #mum");
+ hXA->GetXaxis()->SetTitle("#alpha");
+ hXA->SetStats(0);
+ hXA->FitSlicesY(gs,0,-1,0,"qnr",&harrFXA);
+ TH1* hXA_1 = (TH1*)harrFXA.At(1);
+ hXA_1->SetStats(0);
+ hXA_1->SetDirectory(0);
+ hXA_1->GetYaxis()->SetTitle("#DeltaX_{#mu}(#DeltaZ_{#mu}), #mum");
+ TH1* hXA_2 = (TH1*)harrFXA.At(2);
+ hXA_2->SetStats(0);
+ hXA_2->SetDirectory(0);
+ hXA_2->GetYaxis()->SetTitle("#DeltaX_{#sigma}(#DeltaZ_{#sigma}), #mum");
+
+ //............................................................................................................
+
+ // dX
+
+ Int_t fitStatusX=0;
+
+ TH1* hdx = hXA->ProjectionY("hdx");
+ hdx->GetXaxis()->SetTitle("#DeltaX, #mum");
+ hdx->SetTitle(Form("#DeltaX for pattID %d",pattID));
+ hdx->SetStats(0);
+ pad3->cd();
+ gs->SetParameters(hdx->GetMaximum(),hdx->GetMean(),hdx->GetRMS());
+ if((hdx->GetEntries())<100) fitStatusX = hdx->Fit("gs","ql");
+ else fitStatusX = hdx->Fit("gs","q");
+
+ if(fitStatusX==0){
+ DeltaXmean[pattID]= gs->GetParameter(1);
+ DeltaXmeanErr[pattID]= gs->GetParError(1);
+ DeltaXsigma[pattID]= gs->GetParameter(2);
+ DeltaXsigmaErr[pattID]= gs->GetParError(2);
+ Chi2x[pattID] = gs->GetChisquare();
+ Int_t varNDFx = gs->GetNDF();
+ if(varNDFx>=0)
+ NDFx[pattID] = varNDFx;
+ else{
+ NDFx[pattID]=-1;
+ }
+ }
+ else{
+ DeltaXmean[pattID]=0;
+ DeltaXmeanErr[pattID]=0;
+ DeltaXsigma[pattID]=0;
+ DeltaXsigmaErr[pattID]=0;
+ Chi2x[pattID] = -1;
+ }
+
+ pad3->cd();
+ hdx->Draw();
+
+
+ //.............................................................................................................
+
+
+ // dZ vs alpha
+
+ TObjArray harrFZA;
+ TH2* hZA = (TH2*)harr->At((pattID-LowerLimit)*4+1);
+ hZA->GetYaxis()->SetTitle("#DeltaZ, #mum");
+ hZA->GetXaxis()->SetTitle("#alpha");
+ hZA->SetStats(0);
+ hZA->FitSlicesY(gs,0,-1,0,"qnr",&harrFZA);
+ TH1* hZA_1 = (TH1*)harrFZA.At(1);
+ hZA_1->SetMarkerColor(8);
+ hZA_1->SetLineColor(8);
+ hZA_1->SetStats(0);
+ TH1* hZA_2 = (TH1*)harrFZA.At(2);
+ hZA_2->SetMarkerColor(8);
+ hZA_2->SetLineColor(8);
+ hZA_2->SetStats(0);
+
+ Double_t maxA_1 = TMath::Max(hXA_1->GetMaximum(),hZA_1->GetMaximum());
+ Double_t minA_1 = TMath::Min(hXA_1->GetMinimum(),hZA_1->GetMinimum());
+ Double_t addtorangeA_1 = (maxA_1-minA_1)*0.1;
+ Double_t rangeA_1 = addtorangeA_1*10;
+
+ Double_t maxA_2 = TMath::Max(hXA_2->GetMaximum(),hZA_2->GetMaximum());
+ Double_t minA_2 = TMath::Min(hXA_2->GetMinimum(),hZA_2->GetMinimum());
+ Double_t addtorangeA_2 = (maxA_2-minA_2)*0.1;
+ Double_t rangeA_2 = addtorangeA_2*10;
+
+ hXA_1->SetMaximum(maxA_1+addtorangeA_1);
+ hXA_1->SetMinimum(minA_1-addtorangeA_1);
+
+ hXA_2->SetMaximum(maxA_2+addtorangeA_2);
+ hXA_2->SetMinimum(minA_2-addtorangeA_2);
+
+ pad1->cd();
+ hXA_1->Draw();
+ hZA_1->Draw("same");
+ TLegend* legA_1 = new TLegend(0.8*1.1* TMath::Pi()/2, maxA_1-0.2*rangeA_1, TMath::Pi()/2, maxA_1,"","");
+ legA_1->AddEntry(hXA_1, " #DeltaX_{#mu}", "l");
+ legA_1->AddEntry(hZA_1, " #DeltaZ_{#mu}", "l");
+ legA_1->SetFillColor(0);
+ legA_1->SetBorderSize(0);
+ legA_1->Draw();
+ pad2->cd();
+ hXA_2->Draw();
+ hZA_2->Draw("same");
+ TLegend* legA_2 = new TLegend(0.8*1.1* TMath::Pi()/2, maxA_2-0.2*rangeA_2, TMath::Pi()/2, maxA_2,"","");
+ legA_2->AddEntry(hXA_2, " #DeltaX_{#sigma}", "l");
+ legA_2->AddEntry((TObject*)0, "", "");
+ legA_2->AddEntry(hZA_2, " #DeltaZ_{#sigma}", "l");
+ legA_2->SetFillColor(0);
+ legA_2->SetBorderSize(0);
+ legA_2->Draw();
+
+ //...........................................................................................
+
+
+ // dZ
+
+ Int_t fitStatusZ=0;
+
+ TH1* hdz = hZA->ProjectionY("hdz");
+ hdz->GetXaxis()->SetTitle("#DeltaZ, #mum");
+ hdz->SetTitle(Form("#DeltaZ for pattID %d",pattID));
+ hdz->SetStats(0);
+ pad6->cd();
+ gs->SetParameters(hdz->GetMaximum(),hdz->GetMean(),hdz->GetRMS());
+ if((hdz->GetEntries())<100) fitStatusZ = hdz->Fit("gs","ql");
+ else fitStatusZ = hdz->Fit("gs","q");
+
+ if(fitStatusZ==0){
+ DeltaZmean[pattID]= gs->GetParameter(1);
+ DeltaZmeanErr[pattID]= gs->GetParError(1);
+ DeltaZsigma[pattID]= gs->GetParameter(2);
+ DeltaZsigmaErr[pattID]= gs->GetParError(2);
+ Chi2z[pattID] = gs->GetChisquare();
+ Int_t varNDFz = gs->GetNDF();
+ if(varNDFz>=0)
+ NDFz[pattID] = varNDFz;
+ else{
+ NDFz[pattID]=-1;
+ }
+ }
+ else{
+ DeltaZmean[pattID]=0;
+ DeltaZmeanErr[pattID]=0;
+ DeltaZsigma[pattID]=0;
+ DeltaZsigmaErr[pattID]=0;
+ Chi2z[pattID] = -1;
+ }
+
+ pad6->cd();
+ hdz->Draw();
+
+ //.............................................................................................................
+
+ // dX vs beta
+
+ TObjArray harrFXB;
+ TH2* hXB = (TH2*)harr->At((pattID-LowerLimit)*4+2);
+ hXB->GetYaxis()->SetTitle("#DeltaX, #mum");
+ hXB->GetXaxis()->SetTitle("#beta");
+ hXB->SetStats(0);
+ hXB->FitSlicesY(gs,0,-1,0,"qnr",&harrFXB);
+ TH1* hXB_1 = (TH1*)harrFXB.At(1);
+ hXB_1->SetDirectory(0);
+ hXB_1->SetStats(0);
+ hXB_1->GetYaxis()->SetTitle("#DeltaX_{#mu}(#DeltaZ_{#mu}), #mum");
+ TH1* hXB_2 = (TH1*)harrFXB.At(2);
+ hXB_2->SetDirectory(0);
+ hXB_2->SetStats(0);
+ hXB_2->GetYaxis()->SetTitle("#DeltaX_{#sigma}(#DeltaZ_{#sigma}), #mum");
+
+ //............................................................................................
+
+ //dZ vs beta
+
+ TObjArray harrFZB;
+ TH2* hZB = (TH2*)harr->At((pattID-LowerLimit)*4+3);
+ hZB->GetYaxis()->SetTitle("#DeltaZ, #mum");
+ hZB->GetXaxis()->SetTitle("#alpha");
+ hZB->SetStats(0);
+ hZB->FitSlicesY(gs,0,-1,0,"qnr",&harrFZB);
+ TH1* hZB_1 = (TH1*)harrFZB.At(1);
+ hZB_1->SetMarkerColor(8);
+ hZB_1->SetLineColor(8);
+ hZB_1->SetStats(0);
+ TH1* hZB_2 = (TH1*)harrFZB.At(2);
+ hZB_2->SetMarkerColor(8);
+ hZB_2->SetLineColor(8);
+ hZB_2->SetStats(0);
+
+ Double_t maxB_1 = TMath::Max(hXB_1->GetMaximum(),hZB_1->GetMaximum());
+ Double_t minB_1 = TMath::Min(hXB_1->GetMinimum(),hZB_1->GetMinimum());
+ Double_t addtorangeB_1 = (maxB_1-minB_1)*0.1;
+ Double_t rangeB_1 = 10*addtorangeB_1;
+
+ Double_t maxB_2 = TMath::Max(hXB_2->GetMaximum(),hZB_2->GetMaximum());
+ Double_t minB_2 = TMath::Min(hXB_2->GetMinimum(),hZB_2->GetMinimum());
+ Double_t addtorangeB_2 = (maxB_2-minB_2)*0.1;
+ Double_t rangeB_2 = 10*addtorangeB_2;
+
+
+ hXB_1->SetMaximum(maxB_1+addtorangeB_1);
+ hXB_1->SetMinimum(minB_1-addtorangeB_1);
+
+ hXB_2->SetMaximum(maxB_2+addtorangeB_2);
+ hXB_2->SetMinimum(minB_2-addtorangeB_2);
+
+ pad4->cd();
+ hXB_1->Draw();
+ hZB_1->Draw("same");
+ TLegend* legB_1 = new TLegend(0.8*1.1* TMath::Pi()/2, maxB_1-0.2*rangeB_1, TMath::Pi()/2, maxB_1,"","");
+ legB_1->AddEntry(hXA_1, " #DeltaX_{#mu}", "l");
+ legB_1->AddEntry((TObject*)0, "", "");
+ legB_1->AddEntry(hZA_1, " #DeltaZ_{#mu}", "l");
+ legB_1->SetFillColor(0);
+ legB_1->SetBorderSize(0);
+ legB_1->Draw();
+ pad5->cd();
+ hXB_2->Draw();
+ hZB_2->Draw("same");
+ TLegend* legB_2 = new TLegend(0.8*1.1* TMath::Pi()/2, maxB_2-0.2*rangeB_2, TMath::Pi()/2, maxB_2,"","");
+ legB_2->AddEntry(hXB_2, " #DeltaX_{#sigma}", "l");
+ legB_2->AddEntry((TObject*)0, "", "");
+ legB_2->AddEntry(hZB_2, " #DeltaZ_{#sigma}", "l");
+ legB_2->SetFillColor(0);
+ legB_2->SetBorderSize(0);
+ legB_2->Draw();
+
+ //cnv->Print("patt0.gif");
+
+ //..................................................................................................................
+
+ // Drawing the title
+
+ const char* text = Form("pattID: %d freq: %f nPix: %d nRow: %d nCol: %d xCOG (fraction): %0.2f + (%0.2f) zCOG (fraction): %0.2f + (%0.2f)",
+ pattID,Freq,nPix,nRow,nCol, (*xCentrePix)[pattID],(*xCentreShift)[pattID],(*zCentrePix)[pattID],
+ (*zCentreShift)[pattID]);
+ const char* text1 = Form("#DeltaX_{#mu}: (%.2f #pm %.2f) #mum #DeltaX_{#sigma}: (%.2f #pm %.2f) #mum #DeltaX_{MC-centre}: %.2f #mum #chi^{2}/NDF: %.2f / %d",
+ DeltaXmean[pattID],DeltaXmeanErr[pattID],DeltaXsigma[pattID],DeltaXsigmaErr[pattID], DeltaXmean[pattID]+(*xCentreShift)[pattID]*pitchX, Chi2x[pattID], NDFx[pattID]);
+ const char* text2 = Form("#DeltaZ_{#mu}: (%.2f #pm %.2f) #mum #DeltaZ_{#sigma}: (%.2f #pm %.2f) #mum #DeltaZ_{MC-centre}: %.2f #mum #chi^{2}/NDF: %.2f / %d",
+ DeltaZmean[pattID],DeltaZmeanErr[pattID],DeltaZsigma[pattID],DeltaZsigmaErr[pattID], DeltaZmean[pattID]+(*zCentreShift)[pattID]*pitchZ, Chi2z[pattID], NDFz[pattID]);
+ TLatex title;
+ title.SetTextSize(0.1);
+ title.SetTextAlign(12);
+ titlepad->cd();
+ title.DrawLatex(0.085,0.75,text);
+ title.DrawLatex(0.085,0.5,text1);
+ title.DrawLatex(0.085,0.25,text2);
+
+
+ //Drawing the pattern
+
+ TH2* patternIm = new TH2F("patternIm","", nCol,-0.5,nCol-0.5,nRow,-0.5,nRow-0.5);
+ for (Int_t ir=0; ir<nRow; ir++) {
+ for (Int_t ic=0; ic<nCol; ic++) {
+ if(((TBits*)(*pattDB)[pattID])->TestBitNumber(ir*nCol+ic)){
+ patternIm->Fill(ic,ir);
+ }
+ }
+ }
+
+ patternIm->SetStats(0);
+ padPattern->cd();
+ patternIm->Draw("Acol");
+
+ //Drawing vertical lines on the pattern
+
+ TObjArray vertlinesarray;
+ for(int i=2; i<=nCol; i++){
+ TLine* vline = new TLine(patternIm->GetXaxis()->GetBinLowEdge(i),patternIm->GetYaxis()->GetBinLowEdge(1),
+ patternIm->GetXaxis()->GetBinLowEdge(i),patternIm->GetYaxis()->GetBinLowEdge(nRow+1));
+ vertlinesarray.AddAtAndExpand(vline,i);
+ }
+
+ for(int i=2; i<=nCol; i++){
+ TLine* v1 = (TLine*)vertlinesarray.At(i);
+ padPattern->cd();
+ v1->Draw();
+ }
+
+ //Drawing horizontal lines on the pattern
+
+ TObjArray horlinesarray;
+ for(int i=2; i<=nRow; i++){
+ TLine* hline = new TLine(patternIm->GetXaxis()->GetBinLowEdge(1),patternIm->GetYaxis()->GetBinLowEdge(i),
+ patternIm->GetXaxis()->GetBinLowEdge(nCol+1),patternIm->GetYaxis()->GetBinLowEdge(i));
+ horlinesarray.AddAtAndExpand(hline,i);
+ }
+
+ for(int i=2; i<=nRow; i++){
+ TLine* hor1 = (TLine*)horlinesarray.At(i);
+ padPattern->cd();
+ hor1->Draw();
+ }
+
+ //Drawing MChit and COG
+
+ TMarker* COG = new TMarker((*zCentrePix)[pattID]+(*zCentreShift)[pattID]-0.5,
+ (*xCentrePix)[pattID]+(*xCentreShift)[pattID]-0.5,20);
+ COG->Draw();
+
+ TMarker* MChit = new TMarker((*zCentrePix)[pattID]+(*zCentreShift)[pattID]+(DeltaZmean[pattID])/pitchZ-0.5, //20 is the pixel pitch
+ (*xCentrePix)[pattID]+(*xCentreShift)[pattID]+(DeltaXmean[pattID])/pitchX-0.5,29);
+ MChit->SetMarkerColor(kWhite);
+ MChit->Draw();
+
+ canvas->Print("angles.pdf");
+ canvas->Clear();
+
+ if(pattID<number2dhisto){
+
+ canvasbis->cd();
+ TPad* titlepadbis = new TPad("titlepadbis","",canvas->GetLeftMargin(),0.825,1-canvas->GetRightMargin(),0.95);
+ titlepadbis->SetFillColor(kYellow);
+ titlepadbis->Draw();
+ TPad* pad1bis = new TPad("pad1bis","",canvasbis->GetLeftMargin(),(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin(),0.5,0.8);
+ pad1bis->Draw();
+ TPad* pad2bis = new TPad("pad2bis","",0.5,(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin(),1-canvasbis->GetRightMargin(),0.8);
+ pad2bis->Draw();
+ TPad* pad3bis = new TPad("pad3bis","",canvasbis->GetLeftMargin(),canvasbis->GetBottomMargin(),0.5,(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin());
+ pad3bis->Draw();
+ TPad* pad4bis = new TPad("pad4bis","",0.5,canvasbis->GetBottomMargin(),1-canvasbis->GetRightMargin(),(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin());
+ pad4bis->Draw();
+ TPad* padPatternbis = new TPad("padPatternbis","",0.75,0.825,0.9,0.95);
+ padPatternbis->Draw();
+
+ titlepadbis->cd();
+ title.DrawLatex(0.085,0.75,text);
+ title.DrawLatex(0.085,0.5,text1);
+ title.DrawLatex(0.085,0.25,text2);
+ pad1bis->cd();
+ hXA->Draw("colz");
+ pad2bis->cd();
+ hXB->Draw("colz");
+ pad3bis->cd();
+ hZA->Draw("colz");
+ pad4bis->cd();
+ hZB->Draw("colz");
+ padPatternbis->cd();
+ patternIm->Draw("Acol");
+ for(int i=2; i<=nCol; i++){
+ TLine* v1 = (TLine*)vertlinesarray.At(i);
+ padPatternbis->cd();
+ v1->Draw();
+ }
+ for(int i=2; i<=nRow; i++){
+ TLine* hor1 = (TLine*)horlinesarray.At(i);
+ padPatternbis->cd();
+ hor1->Draw();
+ }
+ COG->Draw();
+ MChit->Draw();
+ canvasbis->Print("angles2D.pdf");
+ canvasbis->Clear();
+ }
+
+ delete patternIm;
+ delete COG;
+ delete MChit;
+
+ horlinesarray.Delete();
+ vertlinesarray.Delete();
+
+}
+
+void UpDateDB(const char* outDBFile){
+
+ printf("\n\n\nStoring data in the DataBase\n\n\n");
+
+ TString outDBFileName = outDBFile;
+ if (outDBFileName.IsNull()) outDBFileName = "clusterTopology2.root";
+
+ TVectorF arrDeltaZmean(nPatterns);
+ TVectorF arrDeltaZmeanErr(nPatterns);
+ TVectorF arrDeltaXmean(nPatterns);
+ TVectorF arrDeltaXmeanErr(nPatterns);
+ TVectorF arrDeltaZsigma(nPatterns);
+ TVectorF arrDeltaZsigmaErr(nPatterns);
+ TVectorF arrDeltaXsigma(nPatterns);
+ TVectorF arrDeltaXsigmaErr(nPatterns);
+ TVectorF arrChi2x(nPatterns);
+ TVectorF arrNDFx(nPatterns);
+ TVectorF arrChi2z(nPatterns);
+ TVectorF arrNDFz(nPatterns);
+
+ for(Int_t ID=0; ID<nPatterns; ID++){
+
+ printf("Processing pattern %d... ", ID);
+
+ arrDeltaZmean[ID]=DeltaZmean[ID];
+ arrDeltaZmeanErr[ID]=DeltaZmeanErr[ID];
+ arrDeltaXmean[ID]=DeltaXmean[ID];
+ arrDeltaXmeanErr[ID]=DeltaXmeanErr[ID];
+ arrDeltaZsigma[ID]=DeltaZsigma[ID];
+ arrDeltaZsigmaErr[ID]=DeltaZsigmaErr[ID];
+ arrDeltaXsigma[ID]=DeltaXsigma[ID];
+ arrDeltaXsigmaErr[ID]=DeltaXsigmaErr[ID];
+ arrChi2x[ID]=Chi2x[ID];
+ arrNDFx[ID]=NDFx[ID];
+ arrChi2z[ID]=Chi2z[ID];
+ arrNDFz[ID]=NDFz[ID];
+
+ printf("done!!\n\n");
+ }
+
+ TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
+ flDB->WriteObject(pattDB,"TopDB","kSingleKey");
+ flDB->WriteObject(NPix,"NPix","kSingleKey");
+ flDB->WriteObject(NRow,"NRow","kSingleKey");
+ flDB->WriteObject(NCol,"NCol","kSingleKey");
+ flDB->WriteObject(pattFR, "TopFreq","kSingleKey");
+ flDB->WriteObject(xCentrePix,"xCOG","kSingleKey");
+ flDB->WriteObject(xCentreShift,"xShift","kSingleKey");
+ flDB->WriteObject(zCentrePix,"zCOG","kSingleKey");
+ flDB->WriteObject(zCentreShift,"zShift","kSingleKey");
+ flDB->WriteObject(&arrDeltaZmean,"DeltaZmean","kSingleKey");
+ flDB->WriteObject(&arrDeltaZmeanErr,"DeltaZmeanErr","kSingleKey");
+ flDB->WriteObject(&arrDeltaXmean,"DeltaXmean","kSingleKey");
+ flDB->WriteObject(&arrDeltaXmeanErr,"DeltaXmeanErr","kSingleKey");
+ flDB->WriteObject(&arrDeltaZsigma,"DeltaZsigma","kSingleKey");
+ flDB->WriteObject(&arrDeltaZsigmaErr,"DeltaZsigmaErr","kSingleKey");
+ flDB->WriteObject(&arrDeltaXsigma,"DeltaXsigma","kSingleKey");
+ flDB->WriteObject(&arrDeltaXsigmaErr,"DeltaXsigmaErr","kSingleKey");
+ flDB->WriteObject(&arrChi2x,"Chi2x","kSingleKey");
+ flDB->WriteObject(&arrChi2z,"Chi2z","kSingleKey");
+ flDB->WriteObject(&arrNDFx,"NDFx","kSingleKey");
+ flDB->WriteObject(&arrNDFz,"NDFz","kSingleKey");
+ //
+ flDB->Close();
+ delete flDB;
+
+ printf("\n\nDB Complete!!\n\n");
+}
\ No newline at end of file
--- /dev/null
+/*
+This macro creates for each pattern in the database (generated by BuildClTop.C) two histograms, for MChit-COG difference in x and z direction respectively
+This histograms are filled with the data from compClusHits.C. With a gaussian fit we get mean value, sigma, chi2, NDF and update the database of patterns
+with this information.
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TArrayI.h"
+#include "TArrayF.h"
+#include "TCanvas.h"
+#include "TBits.h"
+#include "TF1.h"
+#include "../ITS/UPGRADE/AliITSUClusterPix.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSUHit.h"
+#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
+#include "AliITSsegmentation.h"
+#include "AliGeomManager.h"
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "TPaveStats.h"
+
+#endif
+
+TObjArray *pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
+TVectorF* pattFR=0; //frequecy of the pattern in the database
+TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
+TVectorF* zCentrePix=0;
+TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
+TVectorF* zCentreShift=0;
+TVectorF* NPix=0; //Number of fired pixels
+TVectorF* NRow=0; //Number of rows of the pattern
+TVectorF* NCol=0; //Number of columns of the pattern
+
+TArrayF DeltaZmean; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TArrayF DeltaZmeanErr;
+TArrayF DeltaXmean; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TArrayF DeltaXmeanErr;
+TArrayF DeltaZsigma; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TArrayF DeltaZsigmaErr;
+TArrayF DeltaXsigma; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TArrayF DeltaXsigmaErr;
+TArrayF Chi2x; //Chi2 of the gaussian fit over x
+TArrayF Chi2z; //Chi2 of the aussian fut over z
+TArrayI NDFx;
+TArrayI NDFz;
+
+TObjArray histoArr;
+
+/*Int_t err1counter=0;
+Int_t err2counter=0;
+Int_t err3counter=0;*/
+
+enum{kXhisto=0,kZhisto=1};
+
+void LoadDB(const char* fname);
+void CreateHistos(TObjArray* harr);
+TH1* GetHistoX(Int_t pattID, TObjArray* harr);
+TH1* GetHistoZ(Int_t pattID, TObjArray* harr);
+void DrawReport(TObjArray* harr);
+void DrawNP(int pattID, TObjArray* harr);
+void UpDateDB(const char* outDBFile);
+
+typedef struct {
+ Int_t evID;
+ Int_t volID;
+ Int_t lrID;
+ Int_t clID;
+ Int_t nPix;
+ Int_t nX;
+ Int_t nZ;
+ Int_t q;
+ Float_t pt;
+ Float_t eta;
+ Float_t phi;
+ Float_t xyz[3];
+ Float_t dX;
+ Float_t dY;
+ Float_t dZ;
+ Bool_t split;
+ Bool_t prim;
+ Int_t pdg;
+ Int_t ntr;
+ Float_t alpha; // alpha is the angle in y-radius plane in local frame
+ Float_t beta; // beta is the angle in xz plane, taken from z axis, growing counterclockwise
+ Int_t nRowPatt;
+ Int_t nColPatt;
+ Int_t pattID;
+ Float_t freq;
+ Float_t xCen;
+ Float_t zCen;
+ Float_t zShift;
+ Float_t xShift;
+} clSumm;
+
+static clSumm Summ;
+
+Int_t nPatterns=0;
+
+void errClus(){
+
+ //importing data
+
+ LoadDB("clusterTopology.root");
+
+ nPatterns = (pattDB->GetEntriesFast());
+
+ printf("\n\n nPatterns %d\n\n",nPatterns);
+
+ DeltaZmean.Set(nPatterns+100);
+ DeltaZmeanErr.Set(nPatterns+100);
+ DeltaXmean.Set(nPatterns+100);
+ DeltaXmeanErr.Set(nPatterns+100);
+ DeltaZsigma.Set(nPatterns+100);
+ DeltaZsigmaErr.Set(nPatterns+100);
+ DeltaXsigma.Set(nPatterns+100);
+ DeltaXsigmaErr.Set(nPatterns+100);
+ Chi2x.Set(nPatterns+100);
+ Chi2z.Set(nPatterns+100);
+ NDFx.Set(nPatterns+100);
+ NDFz.Set(nPatterns+100);
+
+ TFile *input = new TFile("clInfo.root","read");
+ TTree *clitsu = (TTree*) input->Get("clitsu");
+
+ Int_t ntotCl = (Int_t) clitsu->GetEntries();
+ TBranch *brdZ = clitsu->GetBranch("dZ");
+ brdZ->SetAddress(&Summ.dZ);
+ TBranch *brdX = clitsu->GetBranch("dX");
+ brdX->SetAddress(&Summ.dX);
+ TBranch *brpattID = clitsu->GetBranch("pattID");
+ brpattID->SetAddress(&Summ.pattID);
+
+ CreateHistos(&histoArr);
+
+ for(Int_t iCl=0; iCl<ntotCl; iCl++){
+ clitsu->GetEntry(iCl);
+ GetHistoX(Summ.pattID, &histoArr)->Fill(Summ.dX);
+ GetHistoZ(Summ.pattID, &histoArr)->Fill(Summ.dZ);
+ }
+
+ DrawReport(&histoArr);
+
+ //printf("\n\nNDF is 0 : %d times\n\nNDF and fitstati not corresponding: %d times\n\nNDF not corresponding : %d times\n\n", err3counter,err1counter, err2counter);
+
+ UpDateDB("clusterTopology2.root");
+
+}
+
+void LoadDB(const char* fname){ //load the data base built by BuildClTopDB.C
+
+ // load database
+ TFile* fl = TFile::Open(fname);
+ if(!fl){printf("Could not find %s",fname); exit(1);}
+ pattDB = (TObjArray*)fl->Get("TopDB");
+ pattFR = (TVectorF*)fl->Get("TopFreq");
+ xCentrePix =(TVectorF*)fl->Get("xCOG");
+ zCentrePix =(TVectorF*)fl->Get("zCOG");
+ xCentreShift =(TVectorF*)fl->Get("xShift");
+ zCentreShift =(TVectorF*)fl->Get("zShift");
+ NPix =(TVectorF*)fl->Get("NPix");
+ NCol =(TVectorF*)fl->Get("NCol");
+ NRow =(TVectorF*)fl->Get("NRow");
+}
+
+void CreateHistos(TObjArray* harr){ // creates two histograms for each pattern in the DB
+
+ printf("\n\n Creating histograms... ");
+
+ Int_t histonumber=0;
+
+ if(!harr) harr = &histoArr;
+
+ for(Int_t i=0; i<nPatterns; i++){
+
+ TH1* h = new TH1F(Form("hX%d", i),Form("#DeltaX for pattID %d", i),100,-30,30);
+ h->SetDirectory(0);
+ harr->AddAtAndExpand(h, i);
+ histonumber++;
+ }
+
+ for(Int_t i=0; i<nPatterns; i++){
+
+ TH1* h = new TH1F(Form("hZ%d", i),Form("#DeltaZ for pattID %d", i),100,-30,30);
+ h->SetDirectory(0);
+ harr->AddAtAndExpand(h, i+nPatterns);
+ histonumber++;
+ }
+
+ printf(" %d histograms created, corresponding to %d patterns\n\n", histonumber, histonumber/2);
+}
+
+TH1* GetHistoX(Int_t pattID, TObjArray* harr){
+
+ TH1* h=0;
+ if(!harr) harr = &histoArr;
+
+ h=(TH1*)harr->At(pattID);
+ if (!h) {printf("Unknown histo id=%d\n",pattID); exit(1);}
+ return h;
+}
+
+TH1* GetHistoZ(Int_t pattID, TObjArray* harr){
+
+ Int_t zID=nPatterns+pattID;
+
+ TH1* h=0;
+ if(!harr) harr = &histoArr;
+
+ h=(TH1*)harr->At(zID);
+ if (!h) {printf("Unknown histo id=%d\n",zID); exit(1);}
+ return h;
+}
+
+void DrawReport(TObjArray* harr)
+{
+ // plot all the nPatterns cluster
+ for (int i=0;i<nPatterns;i++) {
+ DrawNP(i,harr);
+ }
+}
+
+void DrawNP(int pattID, TObjArray* harr)
+{
+ static TF1* gs = new TF1("gs","gaus",-500,500);
+ if (!harr) harr = &histoArr;
+
+ printf("\n\nProcessing #DeltaX of pattern %d...",pattID);
+
+ Int_t fitStatusX=0;
+
+ TH1* hdx = (TH1*)harr->At(pattID);
+ hdx->GetXaxis()->SetTitle("#DeltaX, #mum");
+ hdx->SetTitle(Form("#DeltaX for pattID %d",pattID));
+ gs->SetParameters(hdx->GetMaximum(),hdx->GetMean(),hdx->GetRMS());
+ if((hdx->GetEntries())<100) fitStatusX = hdx->Fit("gs","Ql0");
+ else fitStatusX = hdx->Fit("gs","Q0");
+
+ if(fitStatusX==0){
+ Double_t px1 = gs->GetParameter(1);
+ DeltaXmean[pattID]=px1;
+ Double_t ex1 = gs->GetParError(1);
+ DeltaXmeanErr[pattID]=ex1;
+ Double_t px2 = gs->GetParameter(2);
+ DeltaXsigma[pattID]=px2;
+ Double_t ex2 = gs->GetParError(2);
+ DeltaXsigmaErr[pattID]=ex2;
+ Double_t ChiSqx = gs->GetChisquare();
+ Chi2x[pattID] = ChiSqx;
+ Int_t varNDFx = gs->GetNDF();
+ if(varNDFx>=0)
+ NDFx[pattID] = varNDFx;
+ else{
+ /*if(varNDFx==0){
+ printf("\n\nNDF = 0, Chi2 = %e for pattern %d\n\n",ChiSqx,pattID);
+ err3counter++;
+ }*/
+ NDFx[pattID]=-1;
+ }
+ }
+ else{
+ DeltaXmean[pattID]=0;
+ DeltaXmeanErr[pattID]=0;
+ DeltaXsigma[pattID]=0;
+ DeltaXsigmaErr[pattID]=0;
+ Chi2x[pattID] = -1;
+ }
+
+ printf("done!!\n\n");
+
+ printf("Processing #DeltaZ of pattern %d... ",pattID);
+
+ Int_t fitStatusZ=0;
+
+ TH1* hdz = (TH1*)harr->At(pattID+nPatterns);
+ hdz->SetTitle(Form("#DeltaZ for pattID %d",pattID));
+ hdz->GetXaxis()->SetTitle("#DeltaZ, #mum");
+ gs->SetParameters(hdz->GetMaximum(),hdz->GetMean(),hdz->GetRMS());
+ if((hdz->GetEntries())<100) fitStatusZ = hdz->Fit("gs","Ql0");
+ else fitStatusZ = hdz->Fit("gs","Q0");
+
+ if(fitStatusZ==0){
+ Double_t pz1=gs->GetParameter(1);
+ DeltaZmean[pattID]=pz1;
+ Double_t ez1 = gs->GetParError(1);
+ DeltaZmeanErr[pattID]=ez1;
+ Double_t pz2 = gs->GetParameter(2);
+ DeltaZsigma[pattID]=pz2;
+ Double_t ez2 = gs->GetParError(2);
+ DeltaZsigmaErr[pattID]=ez2;
+ Double_t ChiSqz = gs->GetChisquare();
+ Chi2z[pattID]=ChiSqz;
+ Int_t varNDFz = gs->GetNDF();
+ if(varNDFz>=0)
+ NDFz[pattID] = varNDFz;
+ else
+ NDFz[pattID] = -1;
+ }
+ else{
+ DeltaZmean[pattID]=0;
+ DeltaZmeanErr[pattID]=0;
+ DeltaZsigma[pattID]=0;
+ DeltaZsigmaErr[pattID]=0;
+ Chi2z[pattID] = -1;
+ }
+
+ /*if((NDFz[pattID]!=NDFx[pattID])&&fitStatusZ!=fitStatusX){
+ printf("\n\nERR: neither NDF nor FitStatus correspond!\n\n");
+ err1counter++;
+ }
+ else if(NDFz[pattID]!=NDFx[pattID]){
+ printf("\n\nERR: NDFx and NDFz donnot correspond!\n\n");
+ err2counter++;
+ }*/
+
+ printf("done!!\n\n");
+
+}
+
+void UpDateDB(const char* outDBFile){
+
+ printf("\n\n\nStoring data in the DataBase\n\n\n");
+
+ TString outDBFileName = outDBFile;
+ if (outDBFileName.IsNull()) outDBFileName = "clusterTopology2.root";
+
+ TVectorF arrDeltaZmean(nPatterns);
+ TVectorF arrDeltaZmeanErr(nPatterns);
+ TVectorF arrDeltaXmean(nPatterns);
+ TVectorF arrDeltaXmeanErr(nPatterns);
+ TVectorF arrDeltaZsigma(nPatterns);
+ TVectorF arrDeltaZsigmaErr(nPatterns);
+ TVectorF arrDeltaXsigma(nPatterns);
+ TVectorF arrDeltaXsigmaErr(nPatterns);
+ TVectorF arrChi2x(nPatterns);
+ TVectorF arrNDFx(nPatterns);
+ TVectorF arrChi2z(nPatterns);
+ TVectorF arrNDFz(nPatterns);
+
+ for(Int_t ID=0; ID<nPatterns; ID++){
+
+ printf("Processing pattern %d... ", ID);
+
+ arrDeltaZmean[ID]=DeltaZmean[ID];
+ arrDeltaZmeanErr[ID]=DeltaZmeanErr[ID];
+ arrDeltaXmean[ID]=DeltaXmean[ID];
+ arrDeltaXmeanErr[ID]=DeltaXmeanErr[ID];
+ arrDeltaZsigma[ID]=DeltaZsigma[ID];
+ arrDeltaZsigmaErr[ID]=DeltaZsigmaErr[ID];
+ arrDeltaXsigma[ID]=DeltaXsigma[ID];
+ arrDeltaXsigmaErr[ID]=DeltaXsigmaErr[ID];
+ arrChi2x[ID]=Chi2x[ID];
+ arrNDFx[ID]=NDFx[ID];
+ arrChi2z[ID]=Chi2z[ID];
+ arrNDFz[ID]=NDFz[ID];
+
+ printf("done!!\n\n");
+ }
+
+ TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
+ flDB->WriteObject(pattDB,"TopDB","kSingleKey");
+ flDB->WriteObject(NPix,"NPix","kSingleKey");
+ flDB->WriteObject(NRow,"NRow","kSingleKey");
+ flDB->WriteObject(NCol,"NCol","kSingleKey");
+ flDB->WriteObject(pattFR, "TopFreq","kSingleKey");
+ flDB->WriteObject(xCentrePix,"xCOG","kSingleKey");
+ flDB->WriteObject(xCentreShift,"xShift","kSingleKey");
+ flDB->WriteObject(zCentrePix,"zCOG","kSingleKey");
+ flDB->WriteObject(zCentreShift,"zShift","kSingleKey");
+ flDB->WriteObject(&arrDeltaZmean,"DeltaZmean","kSingleKey");
+ flDB->WriteObject(&arrDeltaZmeanErr,"DeltaZmeanErr","kSingleKey");
+ flDB->WriteObject(&arrDeltaXmean,"DeltaXmean","kSingleKey");
+ flDB->WriteObject(&arrDeltaXmeanErr,"DeltaXmeanErr","kSingleKey");
+ flDB->WriteObject(&arrDeltaZsigma,"DeltaZsigma","kSingleKey");
+ flDB->WriteObject(&arrDeltaZsigmaErr,"DeltaZsigmaErr","kSingleKey");
+ flDB->WriteObject(&arrDeltaXsigma,"DeltaXsigma","kSingleKey");
+ flDB->WriteObject(&arrDeltaXsigmaErr,"DeltaXsigmaErr","kSingleKey");
+ flDB->WriteObject(&arrChi2x,"Chi2x","kSingleKey");
+ flDB->WriteObject(&arrChi2z,"Chi2z","kSingleKey");
+ flDB->WriteObject(&arrNDFx,"NDFx","kSingleKey");
+ flDB->WriteObject(&arrNDFz,"NDFz","kSingleKey");
+ //
+ flDB->Close();
+ delete flDB;
+
+ printf("\n\nDB Complete!!\n\n");
+}
--- /dev/null
+/*
+In this macro we create groups of patterns with similar characteristics. First we load the pattern database with information about MC-hit.
+We do not group patterns with a frequency below a treshold (frequencyTreshold).
+We define 7 different ID. 2 for the sigma of the MChit-COG, x and z direction (sigmaXID and sigma ZID);
+2 for the distance between COG and centre of the pixel, x and z (shiftXID and shiftZID);
+2 for the distance between MChit and centre of the pixel, x and z direction (biasXID ad biasZID).
+1 ID for the number of fired pixels.
+We define the binning (number/width of bins) and assign the previous IDs.
+We decide which the grouping method we want to use: sigma, shift and number of pixels (kShift) or sigma, bias and number of pixels (kBias)
+Finally we assign group ID. Frequent patterns (above the treshold), form a one-pattern group. The tohers are in the same group if they have
+the same IDs.
+A fil.txt is print qith he inormation for each apttern, including the IDs and the pattID of the patterns in the same group.
+
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TF1.h"
+#include "TArrayI.h"
+#include "TArrayF.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TPavesText.h"
+#include "TLatex.h"
+#include "TBits.h"
+#include "TGraph.h"
+#include "TStopwatch.h"
+#include "TMath.h"
+#include "../ITS/UPGRADE/AliITSUClusterPix.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSUHit.h"
+#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
+#include "AliITSsegmentation.h"
+#include "AliGeomManager.h"
+
+#endif
+
+TObjArray histoArr;
+TObjArray *pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
+TVectorF* pattFR=0; //frequency of the pattern in the database
+TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
+TVectorF* zCentrePix=0;
+TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
+TVectorF* zCentreShift=0;
+TVectorF* NPix=0;//Number of fired pixels
+TVectorF* NRow=0;//Number of rows of the pattern
+TVectorF* NCol=0;//Number of columns of the pattern
+TVectorF* DeltaZmean=0; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TVectorF* DeltaXmean=0; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TVectorF* DeltaZsigma=0; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TVectorF* DeltaXsigma=0; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TVectorF* DeltaZmeanErr=0;
+TVectorF* DeltaXmeanErr=0;
+TVectorF* DeltaZsigmaErr=0;
+TVectorF* DeltaXsigmaErr=0;
+
+//Defining the ID to create the groups.
+TArrayI sigmaXID;
+TArrayI sigmaZID;
+TArrayI shiftXID;
+TArrayI shiftZID;
+TArrayI biasXID;
+TArrayI biasZID;
+TArrayI NPixID;
+TArrayI groupID;
+Float_t totalGroupFreq=0;
+
+void LoadDB(const char* namefile);
+
+Int_t nPatterns=0;
+
+//Defining the frequency treshold under which ti group patterns
+Float_t frequencyTreshold = 0.001;
+//Defining the bins
+Int_t NumberofSigmaXbins=1000;
+Int_t NumberofSigmaZbins=1000;
+Int_t NumberofShiftXbins=20;
+Int_t NumberofShiftZbins=20;
+Int_t NumberofBiasXbins=2000;
+Int_t NumberofBiasZbins=2000;
+//Defining the boundaries of the bins concerning the number of pixel
+const Int_t limitsnumber = 4;
+//Defining the width of the bins to store patterns with similar sigma/shift/bias (in micron)
+Float_t sigmaXwidth=3;
+Float_t sigmaZwidth=3;
+Float_t biasXwidth=5;
+Float_t biasZwidth=5;
+Float_t shiftXwidth=1./NumberofShiftXbins; //(fraction of the pitch)
+Float_t shiftZwidth=1./NumberofShiftZbins; //(fraction of the pitch)
+
+Int_t nPixlimits[limitsnumber] = {4,10,25,50};
+
+enum{kShift=0, kBias=1};
+//Select kShift to group wrt sigma & COG-centreOfThePixel distance
+//Select kBias to group wrt sigma & MChit-centreOfThePixel distance
+
+Int_t groupingMethod = kShift;
+
+Int_t tempgroupID=0;
+
+void groupdef(){
+
+ //Importing Data
+ LoadDB("clusterTopology2.root");
+
+ //Define moudule segmentation
+
+ AliGeomManager::LoadGeometry("geometry.root");
+ AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+ AliITSUClusterPix::SetGeom(gm);
+ const AliITSsegmentation* segm = gm->GetSegmentation(0);
+ Float_t pitchX = segm->Dpx(0)*10000; // for pitch in X in micron
+ printf("pitchX: %f\n",pitchX);
+ Float_t pitchZ = segm->Dpz(0)*10000; // for pitch in Z in micron
+ printf("pitchX: %f\n",pitchX);
+ delete gm;
+
+ //Setting the number of patterns
+ nPatterns = (pattDB->GetEntriesFast());
+
+ sigmaXID.Set(nPatterns);
+ sigmaZID.Set(nPatterns);
+
+ //Assign -2 to frequent patterns, not to group, and -1 to rare clusters, -3 to patterns with sigma set to zero
+ for(Int_t ID=0; ID<nPatterns; ID++){
+ printf("Assign temporary sigma ID to pattern %d... ",ID);
+ if((*pattFR)[ID]>frequencyTreshold){
+ sigmaXID[ID]=-2;//In order not to consider it within the groups
+ sigmaZID[ID]=-2;//In order not to consider it within the groups
+ }
+ else{
+ sigmaXID[ID]=-1;
+ sigmaZID[ID]=-1;
+ totalGroupFreq+=(*pattFR)[ID];
+ }
+ printf("done\n\n");
+ }
+
+ //Assign to similar patterns the same sigmaXID
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assign sigmaXID to pattern %d... ",i);
+ if(sigmaXID[i]==-1){
+ if((*DeltaXsigma)[i]==0){
+ sigmaXID[i]=-3; // In order not to cosider it within the groups
+ }
+ else{
+ for(Int_t j=0; j<NumberofSigmaXbins; j++){
+ if(j*sigmaXwidth < (*DeltaXsigma)[i] && (*DeltaXsigma)[i]<= (j+1)*sigmaXwidth){
+ sigmaXID[i]=j+1;
+ break;
+ }
+ }
+ }
+ }
+ printf("done!!\n\n");
+ }
+
+ //Assign to similar patterns the same sigmaZID
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assign sigmaZID to pattern %d... ",i);
+ if(sigmaZID[i]==-1){
+ if((*DeltaZsigma)[i]==0){
+ sigmaZID[i]=-3; // In order not to cosider it within the groups
+ }
+ else{
+ for(int j=0; j<NumberofSigmaZbins ; j++){
+ if(j*sigmaZwidth < (*DeltaZsigma)[i] && (*DeltaZsigma)[i]<= (j+1)*sigmaZwidth){
+ sigmaZID[i]=j+1;
+ break;
+ }
+ }
+ }
+ }
+ printf("done!!\n\n");
+ }
+
+ //assigning shiftID
+
+ shiftXID.Set(nPatterns);
+ shiftZID.Set(nPatterns);
+
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assign shiftXID to pattern %d... ",i);
+
+ for(int j=-NumberofShiftXbins/2; j<NumberofShiftXbins/2; j++){
+
+ if(j*shiftXwidth < (*xCentreShift)[i] && (*xCentreShift)[i]<= (j+1)*shiftXwidth){
+ shiftXID[i]=j+1;
+ printf("done!!\n\n");
+ break;
+ }
+ }
+ }
+
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assign shiftZID to pattern %d... ",i);
+
+ for(int j=-NumberofShiftZbins/2; j<NumberofShiftZbins/2; j++){
+ if(j*shiftZwidth < (*zCentreShift)[i] && (*zCentreShift)[i]<= (j+1)*shiftZwidth){
+ shiftZID[i]=j+1;
+ printf("done!!\n\n");
+ break;
+ }
+ }
+ }
+
+ //assigning BiasID
+
+ biasXID.Set(nPatterns);
+ biasZID.Set(nPatterns);
+
+ //Setting all the bias ID to zero
+
+ for(Int_t i=0; i<nPatterns; i++){
+ biasXID[i]=0;
+ biasZID[i]=0;
+ }
+
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assign biasXID to pattern %d... ",i);
+ for(Int_t j=-NumberofBiasXbins/2; j<NumberofBiasXbins/2; j++){
+ if(j*biasXwidth < ((*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX))
+ && ((*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX))<= (j+1)*biasXwidth){
+ biasXID[i]=j+1;
+ break;
+ }
+ }
+ printf("done!!\n\n");
+ }
+
+ for(Int_t i=0; i<nPatterns; i++){
+ biasXID[i]=0;
+ biasZID[i]=0;
+ }
+
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assign biasZID to pattern %d... ",i);
+ for(Int_t j=-NumberofBiasZbins/2; j<NumberofBiasZbins/2; j++){
+ if(j*biasZwidth < ((*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ))
+ && ((*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ))<= (j+1)*biasZwidth){
+ biasZID[i]=j+1;
+ break;
+ }
+ }
+ printf("done!!\n\n");
+ }
+
+
+
+ //Assigning NPixID
+
+ NPixID.Set(nPatterns);
+
+ for(Int_t i=0; i<nPatterns; i++){
+ printf("Assigning NPixID to pattern %d...", i);
+ if((*NPix)[i]<=nPixlimits[0]){
+ NPixID[i]=0;
+ printf("done!!\n\n");
+ }
+ else if(nPixlimits[0]<(*NPix)[i] && (*NPix)[i]<=nPixlimits[limitsnumber-1]){
+ for(Int_t j=0; j<nPatterns; j++ ){
+ if(nPixlimits[j]<(*NPix)[i] && (*NPix)[i]<=nPixlimits[j+1]){
+ NPixID[i]=j+1;
+ printf("done!!\n\n");
+ break;
+ }
+ }
+ }
+ else if((*NPix)[i]>nPixlimits[limitsnumber-1]){
+ NPixID[i]=limitsnumber+1;
+ printf("done!!\n\n");
+ }
+ }
+
+ //Assigning groupID
+
+ groupID.Set(nPatterns);
+
+ //Assign -2 to frequent patterns, not to group, and -1 to rare clusters
+ for(Int_t ID=0; ID<nPatterns; ID++){
+ printf("Assign temporary group ID to pattern %d... ",ID);
+ groupID[ID]=-1;
+ printf("done\n\n");
+ }
+ Int_t k=0;
+ while((*pattFR)[k]>frequencyTreshold){
+ groupID[k]=tempgroupID;
+ tempgroupID++;
+ k++;
+ }
+
+ if(groupingMethod==kShift){
+ for(Int_t i=0; i<nPatterns; i++){
+ if(groupID[i]!=-1) continue;
+ groupID[i]=tempgroupID;
+ printf("Assigning group ID %d... ",tempgroupID);
+ for(Int_t j=i+1; j<nPatterns; j++){
+ if(sigmaXID[j]==-3){
+ groupID[j]=-1;
+ continue;
+ }
+ else if(sigmaXID[j]==sigmaXID[i] && sigmaZID[j]==sigmaZID[i] &&
+ shiftXID[j]==shiftXID[i] && shiftZID[j]==shiftZID[i] &&
+ NPixID[i]==NPixID[j]) groupID[j]=tempgroupID;
+ }
+ printf("done!!\n\n");
+ tempgroupID++;
+ }
+ }
+ else if(groupingMethod==kBias){
+ for(Int_t i=0; i<nPatterns; i++){
+ if(groupID[i]!=-1) continue;
+ groupID[i]=tempgroupID;
+ printf("Assigning group ID %d... ",tempgroupID);
+ for(Int_t j=i+1; j<nPatterns; j++){
+ if(sigmaZID[j]==-3){
+ groupID[j]=-1;
+ continue;
+ }
+ else if(sigmaXID[j]==sigmaXID[i] && sigmaZID[j]==sigmaZID[i]
+ && biasXID[j]==biasXID[i] && biasZID[j]==biasZID[i]
+ && NPixID[i]==NPixID[j]) groupID[j]=tempgroupID;
+ }
+ printf("done!!\n\n");
+ tempgroupID++;
+ }
+ }
+
+ ofstream a("groupdef.txt");
+
+ //setw(55) << "patterns in the group\n" << endl;
+
+ a << Form("A NEGATIVE ID means that the pattern is not in a group.") << endl << endl <<
+ Form("EXCEPTION: biasID CAN be negative") << endl <<
+ "\n\n......................................................................................." <<
+ "................................................................................................\n\n";
+
+ for(int i=0; i<nPatterns; i++){
+
+ printf("Writing info about pattern %d ...", i);
+
+ a << setw(30) << Form("pattID: %d",i) << setw(30) << Form("freq: %f", (*pattFR)[i]) << setw(30) << Form("NPix: %d",Int_t((*NPix)[i]))<< setw(45) << Form("NRow: %d ",Int_t((*NRow)[i])) << setw(45) <<
+ Form("NCol: %d", Int_t((*NCol)[i])) << endl << endl
+ << setw(45) << Form("DeltaXmean: %f (%f)",(*DeltaXmean)[i],(*DeltaXmeanErr)[i]) << setw(45) <<
+ Form("DeltaZmean: %f (%f)",(*DeltaZmean)[i],(*DeltaZmeanErr)[i])<<
+ setw(45) << Form("DeltaXsigma: %f (%f)",(*DeltaXsigma)[i],(*DeltaXsigmaErr)[i]) << setw(45) <<
+ Form("DeltaZsigma: %f (%f)",(*DeltaZsigma)[i],(*DeltaZsigmaErr)[i]) << endl << endl << setw(30) <<
+ Form("xShift: %f",(*xCentreShift)[i]) << setw(45) << Form("zShift: %f",(*zCentreShift)[i]) << setw(45) << Form("BiasX(MChit-centrePix): %f",(*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX)) << setw(45) <<
+ Form("BiasZ(MChit-centrePix): %f",(*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ)) << endl << endl <<
+ setw(30) << Form("sigmaXID: %d",sigmaXID[i]) << setw(15) << Form("sigmaZID: %d",sigmaZID[i]) << setw(15) <<
+ Form("shiftXID: %d",shiftXID[i]) << setw(15) << Form("shiftZID: %d",shiftZID[i]) << setw(15)<< Form("biasXID: %d",biasXID[i]) <<
+ setw(15) << Form("biasZID: %d",biasZID[i]) << setw(15) << Form("NPixID: %d", NPixID[i]) << endl << endl << setw(30) << Form("groupID: %d", groupID[i])
+ << endl << endl << setw(30) << "patterns in this group:\n\n";
+
+ Int_t matchcounter=0;
+
+ for(int j=0; j<nPatterns; j++){
+ if(matchcounter!=0 && (matchcounter%20)==0){
+ a << endl;
+ matchcounter=0;
+ }
+ if(groupID[j]==groupID[i]) {
+ a << j << ", ";
+ matchcounter++;
+ }
+ }
+
+ a << "\n\n......................................................................................." <<
+ "................................................................................................\n\n";
+
+
+ printf("done!\n\n");
+ }
+
+ a.close();
+
+ printf("%d groups found!!!!\n\n",tempgroupID);
+
+ printf("\n\nThe total frequency of the patterns in group is %f\n\n",totalGroupFreq);
+
+ /*
+
+ TVectorF groupDeltaX(tempGID+100);
+ TVectorF groupShiftX(tempGID+100);
+ TVectorF groupDeltaZ(tempGID+100);
+ TVectorF groupShiftZ(tempGID+100);
+ TVectorF patternNum(tempGID+100);
+ groupDeltaX.Zero();
+ groupShiftX.Zero();
+ groupDeltaZ.Zero();
+ groupShiftZ.Zero();
+ patternNum.Zero();
+
+ ofstream b("groupDB.txt");
+
+ b<<setw(15)<<"groupID"<<setw(25)<<"patterns in the group"<<setw(15)<<
+ "DeltaX"<<setw(15)<<"ShiftX"<<setw(15)<<"DeltaZ"<<setw(15)<<"ShiftZ\n"<<endl;
+
+ TCanvas* c = new TCanvas("c","Patterns groups",900,600);
+ TH1F* h = new TH1F("h","Patterns groups",tempGID,-0.5, tempGID-0.5);
+ h->GetXaxis()->SetTitle("group ID");
+ h->SetStats(0);
+
+ for(Int_t gid=0; gid<tempGID; gid++){
+
+ Int_t PatternNumberInGroup=0;
+ Float_t freqSum=0.;//It is the sum of the frequencies of the patterns in the same group
+
+ for(Int_t pattID=0; pattID<nPatterns; pattID++){
+ if (groupID[pattID]==gid)
+ {
+ groupDeltaX[gid]+=(*DeltaXsigma)[pattID]*(*pattFR)[pattID];
+ groupShiftX[gid]+=(*xCentreShift)[pattID]*(*pattFR)[pattID];
+ groupDeltaZ[gid]+=(*DeltaZsigma)[pattID]*(*pattFR)[pattID];
+ groupShiftZ[gid]+=(*zCentreShift)[pattID]*(*pattFR)[pattID];
+ freqSum+=(*pattFR)[pattID];
+ PatternNumberInGroup++;
+
+ h->Fill(gid);
+ }
+ }
+
+ groupDeltaX[gid]=groupDeltaX[gid]/freqSum;
+ groupShiftX[gid]=groupShiftX[gid]/freqSum;
+ groupDeltaZ[gid]=groupDeltaZ[gid]/freqSum;
+ groupShiftZ[gid]=groupShiftZ[gid]/freqSum;
+ patternNum[gid]=PatternNumberInGroup;
+
+ b<<setw(15)<<gid<<setw(25)<<patternNum[gid]<<setw(15)<<
+ groupDeltaX[gid]<<setw(15)<<groupShiftX[gid]<<setw(15)<<
+ groupDeltaZ[gid]<<setw(15)<<groupShiftZ[gid]<<"\n"<<endl;
+ }
+
+ c->cd();
+ h->Draw();
+
+ TPavesText* info = new TPavesText(0.5,1,0.5,1,1,"nb");
+ info->AddText(Form("Number of groups: %d",tempGID));
+ info->AddText("#delta_{X}<5 & #delta_{Z}<5");
+ info->AddText("#DeltaX<0.05 & #DeltaZ<0.05");
+ c->cd();
+ info->Draw();
+
+ b.close();
+ */
+}
+
+void LoadDB(const char* fname){
+
+ printf("\n\nLoading DB... ");
+ // load database
+ TFile* fl = TFile::Open(fname);
+ if(!fl){printf("Could not find %s",fname); exit(1);}
+ pattDB = (TObjArray*)fl->Get("TopDB");
+ pattFR = (TVectorF*)fl->Get("TopFreq");
+ xCentrePix =(TVectorF*)fl->Get("xCOG");
+ zCentrePix =(TVectorF*)fl->Get("zCOG");
+ xCentreShift =(TVectorF*)fl->Get("xShift");
+ zCentreShift =(TVectorF*)fl->Get("zShift");
+ NPix =(TVectorF*)fl->Get("NPix");
+ NCol =(TVectorF*)fl->Get("NCol");
+ NRow =(TVectorF*)fl->Get("NRow");
+ DeltaXmean =(TVectorF*)fl->Get("DeltaXmean");
+ DeltaZmean =(TVectorF*)fl->Get("DeltaZmean");
+ DeltaXsigma =(TVectorF*)fl->Get("DeltaXsigma");
+ DeltaZsigma =(TVectorF*)fl->Get("DeltaZsigma");
+ DeltaXmeanErr =(TVectorF*)fl->Get("DeltaXmeanErr");
+ DeltaZmeanErr =(TVectorF*)fl->Get("DeltaZmeanErr");
+ DeltaZsigmaErr =(TVectorF*)fl->Get("DeltaXsigmaErr");
+ DeltaXsigmaErr =(TVectorF*)fl->Get("DeltaZsigmaErr");
+ printf("done!!\n\n");
+}
\ No newline at end of file
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff