QAmacro_emc_3.C

Go to the documentation of this file.

#include "../auxi.C"
int QAmacro_emc_3()
{
  // Macro loads a file after reconstruction and plots difference between initial direction of particle and angular position of cluster
        // The following part of macro access RunTimeDataBase and initialize PndEmcMapper from it
        FairRunAna *fRun= new FairRunAna();
        fRun->SetInputFile("sim_emc.root");
        fRun->SetOutputFile("dummy_out.root");

        FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
        FairParRootFileIo* parInput1 = new FairParRootFileIo();
        parInput1->open("simparams.root");

        TString emcAsciiPar = gSystem->Getenv("VMCWORKDIR");
        emcAsciiPar += "/macro/params/";
        emcAsciiPar += "emc.par";

        FairParAsciiFileIo* parInput2 = new FairParAsciiFileIo();
        parInput2->open(emcAsciiPar.Data(),"in");

        rtdb->setFirstInput(parInput1);
        rtdb->setSecondInput(parInput2);

        PndEmcGeoPar *geoPar = (PndEmcGeoPar*) rtdb->getContainer("PndEmcGeoPar");
        fRun->Init();

        geoPar->InitEmcMapper();
        std::cout<<""<<std::endl;

        TFile* f = new TFile("full_emc.root"); //file you want to analyse
        TTree *t=(TTree *) f->Get("pndsim") ;

        TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
        t->SetBranchAddress("EmcCluster",&cluster_array);
        TClonesArray* digi_array=new TClonesArray("PndEmcDigi");
        t->SetBranchAddress("EmcDigi",&digi_array);

        TFile* fsim = new TFile("sim_emc.root"); //file you want to analyse
        TTree *tsim=(TTree *) fsim->Get("pndsim") ;

        TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
        tsim->SetBranchAddress("MCTrack",&mctrack_array);

        TVector3 photon_momentum;

        double cluster_energy;
        double cluster_energy_calibrated;
        double cluster_theta, cluster_phi; //position of the cluster
        double theta, phi; // angular position of the initial particle
        double theta_diff, phi_diff;
        int ndigi, npoint;
        double max_energy=0;

    TH1F *ht=new TH1F("ht","Theta distribution",200,0.,180);
        TH1F *h1= new TH1F("h1","Theta difference",200,-5.,5.);
        TH1F *h2= new TH1F("h2","Phi difference",200,-5.,5.);
        TH1F *h3= new TH1F("h3","Cluster energy",100,0.85,1.05);
        TH2F *h2theta= new TH2F("h2theta","Theta difference",200,0.,180.,200,-5.,5.);
        TH2F *h2phi= new TH2F("h2phi","Phi difference",200,0.,180.,200,-5.,5.);
        TH1F *hE1= new TH1F("hE1","E1",200,0.,1.05);
        TH1F *hE1E9= new TH1F("hE1E9","E1 / E9",200,0.,1.05);
        TH1F *hE9E25= new TH1F("hE9E25","E9 / E25",200,0.,1.05);

        // Cluster angular position
        // Entrance point is determined by minimal time

        // Calibrartor version 3 corresponds to non-uniformity of light collection switched on
        PndEmcAbsClusterCalibrator * calibrator1= PndEmcClusterCalibrator::MakeEmcClusterCalibrator(1, 3);

        // Cluster energy
        for (Int_t j=0; j< t->GetEntriesFast(); j++)
        {
                t->GetEntry(j);
                for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
                {
                        PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
                        cluster_energy=cluster->energy();
                        cluster_energy_calibrated=calibrator1->Energy(cluster);
                        if ((cluster->NumberOfDigis()>1)&&(cluster_energy>0.02))
                        h3->Fill(cluster_energy_calibrated);
                        PndEmcClusterEnergySums esum(*cluster, digi_array);
                        hE1->Fill(esum.E1());
                        hE1E9->Fill(esum.E1E9());
                        hE9E25->Fill(esum.E9E25());
                        TVector3 cluster_pos=cluster->where();
                        ht->Fill(cluster_pos.Theta()*180./3.1415);
                }
        }

        for (Int_t j=0; j< t->GetEntriesFast(); j++)//t->GetEntriesFast()
        {
                t->GetEntry(j);
                tsim->GetEntry(j);

                PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
                photon_momentum=mctrack->GetMomentum();
                theta=photon_momentum.Theta();
                phi=photon_momentum.Phi();


                // Loop over clusters
                // If we have 1 initial particle and several cluster
                // we can separate cluster from the first interaction by maximum energy

                max_energy=0;

                for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
                {
                        PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
                        cluster_energy=cluster->energy();
                        if (cluster_energy>max_energy)
                        {
                                max_energy=cluster_energy;
                                TVector3 cluster_pos=cluster->where();
                                cluster_theta=cluster_pos.Theta();
                                cluster_phi=cluster_pos.Phi();
                        }

                }
                if (max_energy>0.6)
                {
                    //cluster_theta-
                    theta_diff=(cluster_theta-theta)*180./TMath::Pi();
                    h1->Fill(theta_diff);
                    h2theta->Fill(theta*TMath::RadToDeg(),theta_diff);
                    //cluster_phi-
                    phi_diff=(cluster_phi-phi)*180./TMath::Pi();
                    h2->Fill(phi_diff);
                    h2phi->Fill(theta*TMath::RadToDeg(),phi_diff);
                }
        }

Bool_t fTest=kTRUE;

Double_t thetaCheckMean=h1->GetMean();
Double_t thetaCheckRMS=h1->GetRMS();

if (TMath::Abs(thetaCheckMean)<0.1 && thetaCheckRMS<0.7 && thetaCheckRMS>0.1)
{
    cout<<"\n Theta Diff - ACCEPTABLE "<<endl;
}
else
{
    cout<<" \n Theta Diff - SOMETHING WENT WRONG "<<endl;
    cout<<"    Theta mean = " << thetaCheckMean << endl;
    cout<<"    Theta RMS  = " << thetaCheckRMS << endl;
    fTest=kFALSE;
}

Double_t phiCheckMean=h2->GetMean();
Double_t phiCheckRMS=h2->GetRMS();

if (TMath::Abs(phiCheckMean)<0.1 && phiCheckRMS<1.0 && phiCheckRMS>0.2)
{
    cout<<"\n Phi Diff - ACCEPTABLE "<<endl;
}
else
{
    cout<<" \n Phi Diff - SOMETHING WENT WRONG "<<endl;
    cout<<"    Phi mean = " << phiCheckMean << endl;
    cout<<"    Phi RMS  = " << phiCheckRMS << endl;
    fTest=kFALSE;
}

Double_t energyCheckMean=h3->GetMean();
Double_t energyCheckRMS=h3->GetRMS();

if (energyCheckMean<1.02 && TMath::Abs(1.0-energyCheckMean)<0.06 && energyCheckRMS<0.04 && energyCheckRMS>0.01)
{
    cout<<"\n Cluster Energy - ACCEPTABLE "<<endl;
}
else
{
    cout<<" \n Cluster Energy - SOMETHING WENT WRONG "<<endl;
    cout<<"    Energy mean = " << energyCheckMean << endl;
    cout<<"    Energy RMS  = " << energyCheckRMS << endl;
    fTest=kFALSE;
}

 if (ht->GetBinContent(15)==0)
   {
     fTest = kFALSE;
     cout << "\n FW endcap absent - BAD" << endl;
   }

if (fTest){
    cout << " Test passed" << endl;
    cout << " All ok " << endl;
}else{
    cout << " Test Failed" << endl;
    cout << " Not Ok " << endl;
}
  CloseGeoManager();
  return 0;
}