// root macro to analyze the clusterization output { int nEvents = 1000; bool verbose = true; // ----- Load libraries ------------------------------------------------ gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C"); basiclibs(); gSystem->Load("libGeoBase"); gSystem->Load("libParBase"); gSystem->Load("libBase"); gSystem->Load("libMCStack"); gSystem->Load("libField"); gSystem->Load("libGen"); gSystem->Load("libPassive"); gSystem->Load("libgenfit"); gSystem->Load("libtrackrep"); gSystem->Load("libtpc"); gSystem->Load("libtpcreco"); gSystem->Load("librecotasks"); gSystem->Load("libMvd"); gSystem->Load("libMvdReco"); // ------------------------------------------------------------------------ // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ // ----- Prepare the histograms ---------------------------------------- TH2D* hisxy = new TH2D("hisxy","",400,-15.,15.,400,-15.,15.); hisxy->SetTitle("MVD MC Hit, xy view;x / cm;y / cm"); TH2D* hisrz = new TH2D("hisrz","",400,-20.,20.,400,-15.,25.); hisrz->SetTitle("MVD MC Hit, rz view;z / cm;r/ cm"); TH2D* hisrphi = new TH2D("hisrphi","",400,-4.,4.,400,-15.,25.); hisrphi->SetTitle("MVD MC Hit, r#phi view;#phi;r/ cm"); TH1D* hisDiffX = new TH1D("hisdiffx","",400,-0.05,0.05); hisDiffX->SetTitle("(MC - RECO) Hit coordinate x;#Deltax / cm;"); TH1D* hisDiffY = new TH1D("hisdiffy","",400,-0.05,0.05); hisDiffY->SetTitle("(MC - RECO) Hit coordinate y;#Deltay / cm;"); TH1D* hisDiffZ = new TH1D("hisdiffz","",400,-0.01,0.01); hisDiffZ->SetTitle("(MC - RECO) Hit coordinate z;#Deltaz / cm;"); TH2D* hisDiffXY = new TH2D("hisdiffxy","",100,-0.05,0.05,100,-0.05,0.05); hisDiffXY->SetTitle("(MC - RECO) Hit coordinates xy view;#Deltax / cm;#Deltay / cm"); TH2D* hisDiffRZ = new TH2D("hisdiffrz","",100,-0.05,0.05,100,-0.00,0.10); hisDiffRZ->SetTitle("(MC - RECO) Hit coordinates rz view;#Deltaz / cm;#Deltar / cm"); TH2D* hisDiffRPhi = new TH2D("hisdiffrphi","",100,-4.,4.,100,-0.00,0.10); hisDiffRPhi->SetTitle("(MC - RECO) Hit coordinates r#phi view;#Delta#phi;#Deltar / cm"); TH1D* hisDiffLocalX = new TH1D("hisdiffLocalx","",400,-0.05,0.05); hisDiffLocalX->SetTitle("(MC - RECO) Local Hit coordinate x;#Deltax_{L} / cm;"); TH1D* hisDiffLocalY = new TH1D("hisdiffLocaly","",400,-0.05,0.05); hisDiffLocalY->SetTitle("(MC - RECO) Local Hit coordinate y;#Deltay_{L} / cm;"); TH1D* hisDiffLocalZ = new TH1D("hisdiffLocalz","",400,-0.001,0.001); hisDiffLocalZ->SetTitle("(MC - RECO) Local Hit coordinate z;#Deltaz_{L} / cm;"); TH2D* hisDiffLocalXY = new TH2D("hisdiffLocalxy","",100,-0.05,0.05,100,-0.05,0.05); hisDiffLocalXY->SetTitle("(MC - RECO) Local Hit XY;#Deltax_{L} / cm;#Deltay_{L} / cm"); TH2D* hisDiffLocalRZ = new TH2D("hisdiffLocalrz","",100,-0.01,0.01,100,-0.00,0.10); hisDiffLocalRZ->SetTitle("(MC - RECO) Local Hit RZ;#Deltaz_{L} / cm;#Deltar_{L} / cm"); TH2D* hisSensorXY = new TH2D("hisLocalxy","",100,-5.,5.,100,-5.,5.); hisSensorXY->SetTitle("Local Sensor Hit XY;x_{L} / cm;y_{L} / cm"); // ------------------------------------------------------------------------ // ----- Get Data from framework -------------------------------------- // load first the file with the GeoManager PndMvdFileNameCreator namecreator("/home/ralfk/mvdmacros/data/lambda0sim.root"); std::string inFile = namecreator.GetSimFileName(false); std::string digiFile = namecreator.GetDigiFileName(false); std::string recoFile = namecreator.GetDigiFileName(false); TFile* f = new TFile(inFile.c_str()); // the sim file you want to analyse TTree *t=(TTree *) f->Get("cbmsim") ; t->AddFriend("cbmsim",digiFile.c_str()); // the digi file t->AddFriend("cbmsim",digiFile.c_str()); // the reco file you want to analyse TClonesArray* mc_array=new TClonesArray("MvdMCPoint"); t->SetBranchAddress("MVDPoint",&mc_array);//Branch names TClonesArray* digiPixel_array=new TClonesArray("MvdDigiPixel"); t->SetBranchAddress("MVDPixelHit",&digiPixel_array);//Branch names TClonesArray* digiStrip_array=new TClonesArray("MvdDigiStrip"); t->SetBranchAddress("MVDStripHit",&digiStrip_array);//Branch names TClonesArray* hit_array=new TClonesArray("MvdHit"); t->SetBranchAddress("MVDHit",&hit_array);//Branch names TGeoManager *geoMan = (TGeoManager*) gDirectory->Get("CBMGeom"); // ------------------------------------------------------------------------ // ----- Start event processing --------------------------------------- TVector3 vecs, vecmc, vecdiff; TVector3 vecslocal, vecmclocal, vecdifflocal; Double_t tmpMaster[3],tmpLocal[3]; TVector2 locals, localmc, localdiff; TGeoHMatrix* currentTransMat; for (Int_t j=0; jGetEntriesFast() ; j++) { t->GetEntry(j); if(mc_array==0){cout<<"No mc_array Nr. "<GetEntriesFast(); i++) { PndMvdHit *hit=(PndMvdHit*)hit_array->At(i); if(hit==0){cout<<"No PndMvdHit Nr. "<GetRefIndex(); // CAUTION this is the point number now! PndMvdMCPoint* point = (PndMvdMCPoint*)mc_array->At(pointid); if(point==0){cout<<"No PndMvdMCPoint Nr. "<cd(point->GetDetName()); vecs.SetXYZ(hit->GetX(), hit->GetY(), hit->GetZ()); vecmc.SetXYZ(0.5*(point->GetX()+point->GetXOut()) ,0.5*(point->GetY()+point->GetYOut()) ,0.5*(point->GetZ()+point->GetZOut()) ); vecdiff = vecmc - vecs; tmpMaster[0]=vecs.x();tmpMaster[1]=vecs.y();tmpMaster[2]=vecs.z(); geoMan->MasterToLocal(tmpMaster,tmpLocal); vecslocal.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]); tmpMaster[0]=vecmc.x();tmpMaster[1]=vecmc.y();tmpMaster[2]=vecmc.z(); geoMan->MasterToLocal(tmpMaster,tmpLocal); vecmclocal.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]); vecdifflocal = vecmclocal - vecslocal; hisxy->Fill(vecs.x(),vecs.y()); // if(vecs.y() > 0.) hisrz->Fill(vecs.z(),vecs.Perp()); // else hisrz->Fill(vecs.z(),-1.*vecs.Perp()); hisrz->Fill(vecs.z(), ((vecs.y()>0.)?1.:-1.)*vecs.Perp() ); hisSensorXY->Fill(vecslocal.x(),vecslocal.y()); hisDiffXY->Fill(vecdiff.x(),vecdiff.y()); hisDiffRZ->Fill(vecdiff.z(),vecdiff.Perp()); hisDiffRPhi->Fill(vecdiff.Phi(),vecdiff.Perp()); hisDiffX->Fill(vecdiff.x()); hisDiffY->Fill(vecdiff.y()); hisDiffZ->Fill(vecdiff.z()); hisDiffLocalX->Fill(vecdifflocal.X()); hisDiffLocalY->Fill(vecdifflocal.Y()); hisDiffLocalZ->Fill(vecdifflocal.Z()); hisDiffLocalXY->Fill(vecdifflocal.X(),vecdifflocal.Y()); hisDiffLocalRZ->Fill(vecdifflocal.Z(),vecdifflocal.Perp()); }//end for i (points in event) }// end for j (events) // ------------------------------------------------------------------------ // ----- Do some fitting ---------------------------------------------- gStyle->SetOptFit(); Double_t par[6]; TF1 *g1 = new TF1("g1","gaus",-0.007,0.007); TF1 *g2 = new TF1("g2","gaus",-0.03,0.03); TF1 *total = new TF1("total","gaus(0)+gaus(3)",-0.03,0.03); total->SetLineColor(kRed); total->SetLineWidth(1); g1->SetLineWidth(1); g2->SetLineWidth(1); hisDiffLocalX->Fit(g1,"R"); hisDiffLocalX->Fit(g2,"R+"); g1->GetParameters(&par[0]); g2->GetParameters(&par[3]); total->SetParameters(par); hisDiffLocalX->Fit(total,"R");//"R+" hisDiffLocalY->Fit(g1,"R"); hisDiffLocalY->Fit(g2,"R+"); g1->GetParameters(&par[0]); g2->GetParameters(&par[3]); total->SetParameters(par); hisDiffLocalY->Fit(total,"R");//"R+" hisDiffLocalZ->Fit(g1,"R"); hisDiffLocalZ->Fit(g2,"R+"); g1->GetParameters(&par[0]); g2->GetParameters(&par[3]); total->SetParameters(par); hisDiffLocalZ->Fit(total,"R");//"R+" /* hisDiffLocalX->Fit("gaus"); hisDiffLocalX->GetFunction("gaus")->SetLineColor(kRed); hisDiffLocalX->GetFunction("gaus")->SetLineWidth(1); hisDiffLocalY->Fit("gaus"); hisDiffLocalY->GetFunction("gaus")->SetLineColor(kRed); hisDiffLocalY->GetFunction("gaus")->SetLineWidth(1); hisDiffLocalZ->Fit("gaus"); hisDiffLocalZ->GetFunction("gaus")->SetLineColor(kRed); hisDiffLocalZ->GetFunction("gaus")->SetLineWidth(1);*/ // ------------------------------------------------------------------------ // ----- Print out into Canvas ---------------------------------------- Int_t a = 3, b = 3; TCanvas* can1 = new TCanvas("MvdTestPlot","MCHit view in MVD",0,0,a*250,b*250); can1->Divide(a,b); can1->cd(1);hisxy->DrawCopy("colz"); can1->cd(2);hisrz->DrawCopy("colz"); can1->cd(3);hisSensorXY->DrawCopy("colz"); can1->cd(4);hisDiffLocalX->DrawCopy(""); can1->cd(5);hisDiffLocalY->DrawCopy(""); can1->cd(6);hisDiffLocalZ->DrawCopy(""); can1->cd(7);hisDiffX->DrawCopy(""); can1->cd(8);hisDiffY->DrawCopy(""); can1->cd(9);hisDiffRPhi->DrawCopy("colz"); can1->Print("outanareco.ps"); // ------------------------------------------------------------------------ // ----- Finish ------------------------------------------------------- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); cout << endl << endl; cout << "Macro finished succesfully." << endl; cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl; cout << endl; // ------------------------------------------------------------------------ }