/* * Trigger.C * * Created on: Feb 12, 2010 * Author: Simone Bianco */ void Trigger(){ // Customize TString MCFile = "Mvd_Test2.root"; TString parFile = "Mvd_Params2.root"; TString recoFile = "Mvd_Test2_reco.root"; const Int_t numScint = 4; // total number of scintillators Int_t thr = 3; // # of scint hits to generate on trigger pulse const Int_t numSens = 6; TString ScintName[numScint]; ScintName[0] = "/cave_1/Vol_0/ScintVol1_0"; ScintName[1] = "/cave_1/Vol_0/ScintVol2_0"; ScintName[2] = "/cave_1/Vol_0/ScintVol3_0"; ScintName[3] = "/cave_1/Vol_0/ScintVol4_0"; TString SensName[numSens]; SensName[0] = "/cave_1/TTVol_0/TTDouble_0/StripActiveTD1_0"; SensName[1] = "/cave_1/TTVol_0/TTDouble_0/StripActiveTD2_0"; SensName[2] = "/cave_1/TTVol_0/TTSingle_0/StripActiveTS3a_0"; SensName[3] = "/cave_1/TTVol_0/TTSingle_0/StripActiveTS3b_0"; SensName[4] = "/cave_1/TTVol_0/TTSingle_0/StripActiveTS4a_0"; SensName[5] = "/cave_1/TTVol_0/TTSingle_0/StripActiveTS4b_0"; // load libs gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C"); gSystem->Load("libTrigger"); // Bool_t Tr[numScint]; Int_t trLevel = 0; Int_t trSucc = 0; // counter for triggers TCanvas *can1 = new TCanvas(); TCanvas *can2 = new TCanvas("hits","hits",50,50,1000,600); can2->Divide(2,2); TH1F *test = new TH1F("Trigger","Trigger",5,-0.5,4.5); TH2F *det1 = new TH2F("det1","det1",100,-0.96,0.96,100,-0.96,0.96); TH2F *det2 = new TH2F("det2","det2",100,-0.96,0.96,100,-0.96,0.96); TH2F *det3 = new TH2F("det3","det3",100,-0.96,0.96,100,-0.96,0.96); TH2F *det4 = new TH2F("det4","det4",100,-0.96,0.96,100,-0.96,0.96); // Load MC info (for the scints) TFile *f = new TFile(MCFile); TFile *prova = new TFile(parFile); prova->Get("FairBaseParSet"); TFile *rec = new TFile(recoFile); TTree *t=(TTree *) f->Get("cbmsim") ; TClonesArray* tr_array=new TClonesArray("PndTriggerMCPoint"); t->SetBranchAddress("TriggerPoint",&tr_array);//Branch names TClonesArray* mc_array=new TClonesArray("PndMCTrack"); t->SetBranchAddress("MCTrack",&mc_array);//Branch names TTree *trec=(TTree *) rec->Get("cbmsim") ; TClonesArray* rec_array=new TClonesArray("PndMvdHit"); trec->SetBranchAddress("MVDHitsStrip",&rec_array);//Branch names if ((t->GetEntries())!=(trec->GetEntries())) cout << "FILES NOT COMPATIBLE!!!" << endl; cout << "MC" << t->GetEntries() << endl; cout << "reco" << trec->GetEntries() << endl; Double_t buff3x,buff3y,buff4x,buff4y; Bool_t fill3, fill4; fGeoH = new PndMvdGeoHandling(gGeoManager); for (Int_t j = 0 ; j < t->GetEntries() ; j++) { t->GetEvent(j); if (j%100 == 0) cout << "Ev. " << j << endl; trec->GetEvent(j); trLevel = 0; fill3 = kFALSE; fill4 = kFALSE; for (Int_t a = 0 ; a < numScint ; a++) { Tr[a]= kFALSE; } for (Int_t y = 0 ; y < tr_array->GetEntries() ; y++) { PndTriggerMCPoint *point = (PndTriggerMCPoint*)tr_array->At(y); PndMCTrack *track = (PndMCTrack*)mc_array->At(point->GetTrackID()); //cout << "Event " << j << ", " << fGeoH->GetPath(point->GetDetName()) << ", point " << y << ", pdg: " << track->GetPdgCode() << endl; if (!((track->GetMotherID())<=0)) { //cout << "Skipping!" << endl; continue; } for (Int_t w = 0 ; w < numScint ; w++) { if ((fGeoH->GetPath(point->GetDetName())) == ScintName[w]) Tr[w] = kTRUE; } } // end loop on hit regarding one event for (Int_t z = 0 ; z < numScint ; z++) { if (Tr[z]) trLevel++; } //cout << "Trigger: " << trLevel << ":" << Tr[0] << "|" << Tr[1] << "|" << Tr[2] << "|" << Tr[3] << "|" << endl; test->Fill(trLevel); // check on trigger condition if (trLevel >= thr) { trSucc++; for(Int_t lrec = 0 ; lrec < rec_array->GetEntries() ; lrec++) { // load hits from clustering PndMvdHit *hit = (PndMvdHit*) rec_array->At(lrec); //cout << fGeoH->GetPath(hit->GetDetName()) << endl; if ( (fGeoH->GetPath(hit->GetDetName())) == SensName[0]) { //cout << "x:" << hit->GetX() << "y" << hit->GetY() << endl; det1->Fill(hit->GetX(),hit->GetY()); } if ( (fGeoH->GetPath(hit->GetDetName())) == SensName[1]) { //cout << "x:" << hit->GetX() << "y" << hit->GetY() << endl; det2->Fill(hit->GetX(),hit->GetY()); } if ( (fGeoH->GetPath(hit->GetDetName())) == SensName[2]) { buff3x = hit->GetX(); if (fill3) { det3->Fill(buff3x,buff3y); //cout << "x:" << hit->GetX() << "y" << hit->GetY() << endl; } else fill3 = kTRUE; } if ( (fGeoH->GetPath(hit->GetDetName())) == SensName[3]) { buff3y = hit->GetY(); if (fill3) { det3->Fill(buff3x,buff3y); //cout << "x:" << hit->GetX() << "y" << hit->GetY() << endl; } else fill3 = kTRUE; } if ( (fGeoH->GetPath(hit->GetDetName())) == SensName[4]) { buff4x = hit->GetX(); if (fill4) { det4->Fill(buff4x,buff4y); //cout << "x:" << hit->GetX() << "y" << hit->GetY() << endl; } else fill4 = kTRUE; } if ( (fGeoH->GetPath(hit->GetDetName())) == SensName[5]) { buff4y = hit->GetY(); if (fill4) { det4->Fill(buff4x,buff4y); //cout << "x:" << hit->GetX() << "y" << hit->GetY() << endl; } else fill4 = kTRUE; } } } } // end loop on events can1->cd(); test->Draw(); Int_t m1,m2,m3,m4,max; m1 = det1->GetBinContent(det1->GetMaximumBin()); m2 = det2->GetBinContent(det2->GetMaximumBin()); m3 = det3->GetBinContent(det3->GetMaximumBin()); m4 = det4->GetBinContent(det4->GetMaximumBin()); max = maximum(m1,m2,m3,m4); can2->cd(1); det1->SetMaximum(max); det1->Draw("COLZ"); can2->cd(2); det2->SetMaximum(max); det2->Draw("COLZ"); can2->cd(3); det3->SetMaximum(max); det3->Draw("COLZ"); can2->cd(4); det4->SetMaximum(max); det4->Draw("COLZ"); cout << "Number of trigger signals: " << trSucc << endl; } Int_t maximum( Int_t x, Int_t y, Int_t z , Int_t w) { Int_t max = x; if ( y > max ) max = y; if ( z > max ) max = z; if ( w > max ) max = w; return max; } //end function maximum //Int_t vols; // //(gGeoManager->GetListOfVolumes())->Print(); // //vols = gGeoManager->GetNNodes(); //cout << vols << endl; // // // // //for (Int_t w = 0 ; w GetVolume(w); // if (!vol->IsAssembly()) // { // char *nome = vol->GetName(); //// cout << vol->GetName() << endl; // cout << nome << endl; // cout << vol->GetUniqueID() << endl; //// gGeoManager->CdNode(w); //// cout << gGeoManager->GetPath() << endl; // } //// cout << (gGeoManager->GetVolume(w))->GetName() << endl; //} // // // //}