// // Calorimeter useful functions and subroutines - Emiliano Mocchiutti // // CaloFunctions.h version 4.00 (2005-11-29) // // Programs in this file are called by other programs and cannot be run by hand. // // Changelog: // // 3.10 - 4.00 (2005-11-29): preparing for the final release, changed checkemilib (PAM_YODALIB -> PAM_YODA). // // 3.09 - 3.10 (2005-11-15): changed langaupro since it did not work in the compiled version; changed WhatToDo to work in compiled version. // // 3.08 - 3.09 (2005-11-11): small changes to PrintFigure subroutine // // 3.07 - 3.08 (2005-10-13): classes are now in another file (caloclasses.h) // // 3.06 - 3.07 (2005-10-10): added subroutine ColorTOFMIP for TOF colors. // // 3.05 - 3.06 (2005-10-04): added subroutine getfileLEVname. // // 3.04 - 3.05 (2005-08-04): Fixed bugs when running on a 64 bit arch. // // 3.03 - 3.04 (2005-07-14): Bug in getLEVname, fixed. // // 3.02 - 3.03 (2005-07-13): Finally CaloFunctions.h can be compiled under ROOT! Added clones of routines in the yodautility file no needed anymore. // // 3.01 - 3.02 (2005-07-11): Some more tuning for the macro compilation. // // 3.00 - 3.01 (2005-07-06): changed getLEVname since it gave some problems when called more than once in a ROOT sesssion. Tracker structurs now placed in a tracker // include file. // // 2.27 - 3.00 (2005-07-05): Started the tuning for the compilation of the macros. Changed order of functions, added include files for ROOT routines, fixed some // implicit casting. Fixed redefinition of class constructors. Small change in the calib structure. // // 2.26 - 2.27 (2005-06-23): Changed class CalorimeterLevel2 and structures to add good, perr,swerr and crc infos. // // 2.25 - 2.26 (2005-06-16): Changed class CalorimeterLevel2 and structures to convert OBT, headc and evfile to integers. // // 2.24 - 2.25 (2005-06-07): Added subroutine checkifpretyodaversion to check if file is older than 050515_007 (since that file there is a new format in TMTC data). // // 2.23 - 2.24 (2005-06-06): Small changes in the struct "calib". Change type of CaloPede subroutine, now it is an integer and not void: return 3 if it doesn't find the // calibration file. // // 2.22 - 2.23 (2005-05-30): Added checkmerging and fillmerging subroutines to avoid processing twice the same event. // // 2.21 - 2.22 (2005-05-19): Speed-up and better behaviour of CaloFindcalibs. No need to change calls to this subroutine, all internal changes. // // 2.20 - 2.21 (2005-05-17): Changed CaloFindCalibs in order to check also the calibration file provided if any but first check the current file. // // 2.19 - 2.20 (2005-05-06): Added mip to structure calib. // // 2.18 - 2.19 (2005-05-05): Added calorimeter class CalorimeterLevel2 and structures CaLevel1 CaLevel2. Changed getLEVname subroutines (still compatible with older version). // // 2.17 - 2.18 (2005-04-19): Added the old subroutine Calo1stcalib with the name Calo1stcalibOLD (to be used in ShowEvent.c). // // 2.16 - 2.17 (2005-04-14): Added classes for calorimeter ADC to MIP conversion. Bug fixed in CalorimeterCalibration class. CaloFindBaseRawNC subroutine added to // calculate baseline from raw data without any compression. Added getEmiFile subroutine to find YODA files (equivalent of getFile in // yodautility.c); sintax is similar to getFile but it doesn't crash after 1111 or 2222 files. // // 2.15 - 2.16 (2005-04-11): Bugs fixes in the CaloFindCalibs, fetchpreviousfile and whatnamewith subroutines. Consider the case of an OBT jump in data. // // 2.14 - 2.15 (2005-04-07): Changes in the CaloFindCalibs subroutine. Added CaloCompressData subroutine, some small bugs fixed. // // 2.13 - 2.14 (2005-04-06): Changes in the CaloFindCalibs subroutine. Added fetchpreviousfile and whatnamewith subroutines. Now // the program looks for calorimeter calibrations also in previous files and associates them to data in // a correct way. // // 2.12 - 2.13 (2005-03-17): Changes in the WhatToDo subroutine. // // 2.11 - 2.12 (2005-03-09): PrintFigure subroutine changed to accomplish the new "saveas" script style. // // 2.10 - 2.11 (2005-03-04): Changed structure CalorimeterCalibration to save also fit parameters. // // 2.09 - 2.10 (2005-03-02): Show calibration CRC errors (in CaloFindCalibs). Added langaupro subroutines. // // 2.08 - 2.09 (2005-02-28): Calorimeter calibration classes changed (required by the new CaloADC2MIP). // // 2.07 - 2.08 (2005-02-25): check the existence of calibration header file in CaloFindCalibs! // // 2.06 - 2.07 (2005-02-24): bug in CaloFindBase! DO NOT put to zero the ADC values if you can't determine the baseline! // // 2.05 - 2.06 (2005-02-24): CalorimeterLevel1 class must be cleared during initialization to avoid the so-called "gundam-bug". // // 2.04 - 2.05 (2005-02-24): Changed variable definition from C/C++ style to ROOT style (int->Int_t). // // 2.03 - 2.04 (2005-02-23): Mispelling in the WhatToDo subroutine, fixed. // // 2.02 - 2.03 (2005-02-18): Calo1stCalib call CaloPede also if we don't have all four the calorimeter calibrations. // // 2.01 - 2.02 (2005-02-10): "progressive number" instead of "event number" in WhatToDo subroutine. // // 2.00 - 2.01 (2005-02-07): Added langau fit subroutines and the PrintFigure subroutine (used in CaloPLANES mainly). Some bug fixed, // added user interaction routine and getfilename and getlevname subroutines to retrieve file name and to give // the level name correctly. // // 1.00 - 2.00 (2005-01-17): Cleanup of the code // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include typedef struct Chmerge { Int_t fromflag; Int_t scanflag; Int_t pscuOBT; Int_t pscucount; Int_t lastOBT[1000]; Int_t lastcount[1000]; } chmerge; typedef struct Evento { Int_t iev; Int_t stwerr[4]; Float_t emin; Float_t perror[4]; Float_t dexy[2][22][96]; Float_t dexyc[2][22][96]; Double_t base[2][22][6]; Float_t calselftrig[4][7]; Float_t calIItrig[4]; Float_t calstriphit[4]; Float_t calDSPtaberr[4]; Float_t calevnum[4]; Float_t shift; } evento; typedef struct Calib { TString tablename; TString basepath; TString yodalev; char *db; Int_t DW; Int_t mysql; Int_t obtjump; Int_t time[4][51]; Int_t ttime[4][51]; Int_t fcode[4][51]; Int_t iev; Int_t cstwerr[4]; Float_t ispaw; Float_t cperror[4]; Float_t mip[2][22][96]; Float_t calped[2][22][96]; Float_t calgood[2][22][96]; Float_t calthr[2][22][6]; Float_t calrms[2][22][96]; Float_t calbase[2][22][6]; Float_t calvar[2][22][6]; Float_t calpuls[2][22][96]; Double_t sbase[2][22][6]; Double_t al_p[5][2]; Int_t trkchi2; Bool_t good2; } calib; void emicheckLib(){ if (!TClassTable::GetDict("pamela::TmtcEvent")){ char *pamyodalib = gSystem->ExpandPathName("$PAM_YODA"); stringstream libYoda; libYoda.str(""); libYoda << pamyodalib << "/lib/libyoda.so"; gSystem->Load(libYoda.str().c_str()); } } void emidigForFiles(TList& out, TSystemDirectory *tsd, string defin){ const char *initialDir = gSystem->pwd(); TList *lnk = tsd->GetListOfFiles(); if (lnk==0) return; TSystemFile *file = (TSystemFile*)lnk->First(); if (file->IsZombie()) { lnk->Delete(); gSystem->cd(initialDir); return; } string *fileDes = new string(file->GetName()); unsigned int loc = 0; while(file){ fileDes = new string(file->GetName()); if(!((fileDes->compare(".") == 0) || (fileDes->compare("..") == 0))){ if (file->IsDirectory()){ emidigForFiles(out, (TSystemDirectory*)file, defin); } else { loc = fileDes->find(defin.c_str(), 0); if (loc != string::npos){ out.AddLast((TObject*)file); } } } file = (TSystemFile*)lnk->After(file); } gSystem->cd(initialDir); } void emimakeAllFriend(TTree* out, TList* input){ TList *keys; TKey *key; TTree *tr; string *base; string *className; unsigned int loc; bool cont = false; TSystemFile *tsf = (TSystemFile*)input->First(); if (tsf->IsZombie()) return; while(tsf){ base = new string(tsf->GetTitle()); base->append("/"); base->append(tsf->GetName()); TFile *tf = new TFile(base->c_str()); if (tf->IsZombie()) { tsf = (TSystemFile*)input->After(tsf); continue; } keys = tf->GetListOfKeys(); key = (TKey*)keys->First(); while(key){ className = new string(key->GetClassName()); loc = className->find("Tree", 0); if(loc != string::npos){ tr = (TTree*)key->ReadObj(); if (tr->GetEntries() > 0){ if(cont){ tr = (TTree*)key->ReadObj(); out->AddFriend(tr->GetName(), tf); } else { out = ((TTree*)key->ReadObj()); cont = true; } } } key = (TKey*)keys->After(key); } tsf = (TSystemFile*)input->After(tsf); } } void emimakeAllChained(TChain& out, TList* input){ TList *keys; TKey *key; TTree *tr; string *base; string *className; unsigned int loc; TSystemFile *tsf = (TSystemFile*)input->First(); if (tsf->IsZombie()) return; while(tsf){ base = new string(tsf->GetTitle()); base->append("/"); base->append(tsf->GetName()); TFile *tf = new TFile(base->c_str()); if (tf->IsZombie()) { tsf = (TSystemFile*)input->After(tsf); continue; } keys = tf->GetListOfKeys(); key = (TKey*)keys->First(); while(key){ className = new string(key->GetClassName()); loc = className->find("Tree", 0); if(loc != string::npos){ tr = (TTree*)key->ReadObj(); if (tr->GetEntries() > 0){ out.Add(base->c_str()); } } key = (TKey*)keys->After(key); } tsf = (TSystemFile*)input->After(tsf); } } TFile* emigetFile(TString base, TString packetType, TString subType = "Event"){ TSystemDirectory *targetDir = new TSystemDirectory("", base); TList *filesList = new TList; TSystemFile *tsf = 0; string *tmpString = 0; tmpString = new string("."); tmpString->append(packetType); tmpString->append("."); tmpString->append(subType); tmpString->append("."); emidigForFiles(*filesList, targetDir, tmpString->c_str()); if ( filesList->IsEmpty() ) { filesList->Delete(); return 0; }; targetDir = new TSystemDirectory("", base); tsf = (TSystemFile*)filesList->First(); tmpString = new string(tsf->GetTitle()); tmpString->append("/"); tmpString->append(tsf->GetName()); targetDir->Delete(); filesList->Delete(); return new TFile(tmpString->c_str()); } Int_t emigetLastNotZeroBin(TH1 *histo){ Int_t lastBin = 0; Int_t tempBin = 0; Stat_t minVal = 0; Stat_t tempVal = 0; Int_t range = histo->GetNbinsX(); tempBin = histo->GetMaximumBin(); tempVal = histo->GetBinContent(tempBin); minVal = tempVal; lastBin = tempBin; for (Int_t i = histo->GetMaximumBin() + 1; i < range; i++){ tempVal = histo->GetBinContent(i); if ((tempVal != 0) && (tempVal < minVal)) { minVal = tempVal; lastBin = i; } } return (Int_t)(lastBin*1.10); } Int_t emigetFirstNotZeroBin(TH1 *histo){ Int_t lastBin = 0; Int_t tempBin = 0; Stat_t minVal = 0; Stat_t tempVal = 0; tempBin = histo->GetMaximumBin(); tempVal = histo->GetBinContent(tempBin); minVal = tempVal; lastBin = tempBin; for (Int_t i = histo->GetMaximumBin() - 1; i > 0; i--){ tempVal = histo->GetBinContent(i); if ((tempVal != 0) && (tempVal < minVal)) { minVal = tempVal; lastBin = i; } } return (Int_t)(lastBin*0.90); } void stringcopy(TString& s1, const TString& s2, Int_t from=0, Int_t to=0){ if ( to == 0 ){ Int_t t2length = s2.Length(); s1 = ""; to = t2length; }; for (Int_t i = from; iIsFileInIncludePath(newfile) ){ return 0; } else { return new TFile(newfile); }; } void CaloCompressData(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ Float_t minstrip = 100000.; Float_t rms = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) { minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e]; rms = calib.calthr[l][m][pre]; }; }; // // compression // if ( minstrip < evento.base[l][m][pre] && minstrip != 100000.){ for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.; }; }; } void CaloFindBase(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ Float_t ominstrip = 100000.; Float_t orms = 0.; Float_t minstrip = 100000.; Float_t rms = 0.; evento.base[l][m][pre] = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) { ominstrip = minstrip; minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e]; orms = rms; rms = calib.calthr[l][m][pre]; }; if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < ominstrip && evento.dexyc[l][m][e]-calib.calped[l][m][e] >minstrip ) { ominstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e]; orms = calib.calthr[l][m][pre]; }; }; if ( minstrip != 100000. ) { if ( ominstrip != 10000.) { minstrip = ominstrip; rms = orms; }; Float_t strip6s = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( (evento.dexyc[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexyc[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) { strip6s += 1.; evento.base[l][m][pre] += (evento.dexyc[l][m][e] - calib.calped[l][m][e]); }; // // compression // if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.; }; if ( strip6s >= 9. ){ Double_t arro = evento.base[l][m][pre]/strip6s; Float_t deci = 1000.*((float)arro - float(int(arro))); if ( deci < 500. ) { arro = double(int(arro)); } else { arro = 1. + double(int(arro)); }; evento.base[l][m][pre] = arro; } else { evento.base[l][m][pre] = 31000.; }; } else { evento.base[l][m][pre] = 31000.; }; } void CaloFindBaseRaw(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ Float_t minstrip = 100000.; Float_t rms = 0.; evento.base[l][m][pre] = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) { minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e]; rms = calib.calthr[l][m][pre]; }; }; if ( minstrip != 100000. ) { Float_t strip6s = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) { strip6s += 1.; evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]); }; // // compression // if ( abs((int)(evento.dexy[l][m][e]-calib.calped[l][m][e])) <= (minstrip+rms) ) { evento.dexyc[l][m][e] = 0.; } else { evento.dexyc[l][m][e] = evento.dexy[l][m][e]; }; }; if ( strip6s >= 9. ){ Double_t arro = evento.base[l][m][pre]/strip6s; Float_t deci = 1000.*((float)arro - float(int(arro))); if ( deci < 500. ) { arro = double(int(arro)); } else { arro = 1. + double(int(arro)); }; evento.base[l][m][pre] = arro; } else { evento.base[l][m][pre] = 31000.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ evento.dexyc[l][m][e] = evento.dexy[l][m][e]; }; }; } else { evento.base[l][m][pre] = 31000.; }; } void CaloFindBaseRawNC(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){ Float_t minstrip = 100000.; Float_t rms = 0.; evento.base[l][m][pre] = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) { minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e]; rms = calib.calthr[l][m][pre]; }; }; if ( minstrip != 100000. ) { Float_t strip6s = 0.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) { strip6s += 1.; evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]); }; // // compression // evento.dexyc[l][m][e] = evento.dexy[l][m][e]; }; if ( strip6s >= 9. ){ Double_t arro = evento.base[l][m][pre]/strip6s; Float_t deci = 1000.*((float)arro - float(int(arro))); if ( deci < 500. ) { arro = double(int(arro)); } else { arro = 1. + double(int(arro)); }; evento.base[l][m][pre] = arro; } else { evento.base[l][m][pre] = 31000.; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ evento.dexyc[l][m][e] = evento.dexy[l][m][e]; }; }; } else { evento.base[l][m][pre] = 31000.; }; } int CaloPede(TString filename, Int_t s, Int_t atime, Calib & calib){ TFile *calheadFile = 0; calheadFile = getEmiFile(filename, "CalibCalPed", "Header"); TFile *calcalibFile = 0; calcalibFile = getEmiFile(filename, "CalibCalPed","Event"); // const char *filen = 0; filen = filename; if ( !calheadFile || !calcalibFile ) { if ( calheadFile ) calheadFile->Close(); if ( calcalibFile ) calcalibFile->Close(); printf("\n\nERROR: not able to open file: %s \n",filen); return(3); }; TTree *tr = (TTree*)calheadFile->Get("Pscu"); tr->AddFriend("CalibCalPed", calcalibFile); pamela::CalibCalPedEvent *ce = 0; pamela::PscuHeader *cph = 0; pamela::EventHeader *ceh = 0; tr->SetBranchAddress("Header", &ceh); tr->SetBranchAddress("Event", &ce); Long64_t ncalibs = tr->GetEntries(); for (Int_t ci = 0; ci < ncalibs ; ci++){ tr->GetEntry(ci); cph = ceh->GetPscuHeader(); if ( atime == cph->GetOrbitalTime()){ calib.iev = ce->iev; if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { for ( Int_t d=0 ; d<11 ;d++ ){ Int_t pre = -1; for ( Int_t j=0; j<96 ;j++){ if ( j%16 == 0 ) pre++; if ( s == 2 ){ calib.calped[0][2*d+1][j] = ce->calped[3][d][j]; calib.cstwerr[3] = ce->cstwerr[3]; calib.cperror[3] = ce->cperror[3]; calib.calgood[0][2*d+1][j] = ce->calgood[3][d][j]; calib.calthr[0][2*d+1][pre] = ce->calthr[3][d][pre]; calib.calrms[0][2*d+1][j] = ce->calrms[3][d][j]; calib.calbase[0][2*d+1][pre] = ce->calbase[3][d][pre]; calib.calvar[0][2*d+1][pre] = ce->calvar[3][d][pre]; }; if ( s == 3 ){ calib.calped[0][2*d][j] = ce->calped[1][d][j]; calib.cstwerr[1] = ce->cstwerr[1]; calib.cperror[1] = ce->cperror[1]; calib.calgood[0][2*d][j] = ce->calgood[1][d][j]; calib.calthr[0][2*d][pre] = ce->calthr[1][d][pre]; calib.calrms[0][2*d][j] = ce->calrms[1][d][j]; calib.calbase[0][2*d][pre] = ce->calbase[1][d][pre]; calib.calvar[0][2*d][pre] = ce->calvar[1][d][pre]; }; if ( s == 0 ){ calib.calped[1][2*d][j] = ce->calped[0][d][j]; calib.cstwerr[0] = ce->cstwerr[0]; calib.cperror[0] = ce->cperror[0]; calib.calgood[1][2*d][j] = ce->calgood[0][d][j]; calib.calthr[1][2*d][pre] = ce->calthr[0][d][pre]; calib.calrms[1][2*d][j] = ce->calrms[0][d][j]; calib.calbase[1][2*d][pre] = ce->calbase[0][d][pre]; calib.calvar[1][2*d][pre] = ce->calvar[0][d][pre]; }; if ( s == 1 ){ calib.calped[1][2*d+1][j] = ce->calped[2][d][j]; calib.cstwerr[2] = ce->cstwerr[2]; calib.cperror[2] = ce->cperror[2]; calib.calgood[1][2*d+1][j] = ce->calgood[2][d][j]; calib.calthr[1][2*d+1][pre] = ce->calthr[2][d][pre]; calib.calrms[1][2*d+1][j] = ce->calrms[2][d][j]; calib.calbase[1][2*d+1][pre] = ce->calbase[2][d][pre]; calib.calvar[1][2*d+1][pre] = ce->calvar[2][d][pre]; }; }; }; }; }; }; tr->Delete(); if ( calheadFile ) calheadFile->Close(); if ( calcalibFile ) calcalibFile->Close(); return(0); } TString getFilename(const TString filename){ const string fil = (const char*)filename; Int_t posiz = fil.find("dw_"); if ( posiz == -1 ) posiz = fil.find("DW_"); if ( posiz == -1 ) return 0; Int_t posiz2 = posiz+13; TString file2; stringcopy(file2,filename,posiz,posiz2); TString pdat(".dat"); stringappend(file2,pdat); return file2; } Int_t fetchpreviousfile(TString ffile, Int_t s){ const char *file = ffile; stringstream fil; fil.str(""); fil << file ; Int_t posiz = fil.str().find("dw_"); Int_t upper = 0; if ( posiz == -1 ) { posiz = fil.str().find("DW_"); upper = 1; }; if ( posiz == -1 ) return(1); // TString nomefile = getFilename(ffile); const char *ufnome = nomefile; stringstream nyear; stringstream nmonth; stringstream nday; stringstream nfno; TString year; TString month; TString day; TString fno; stringcopy(year, ufnome, 3, 5); stringcopy(month, ufnome, 5, 7); stringcopy(day, ufnome, 7, 9); stringcopy(fno, ufnome, 10, 13); const char *cyear = year; const char *cmonth = month; const char *cday = day; const char *cfno = fno; Int_t jy = atoi(cyear); Int_t jm = atoi(cmonth); Int_t jd = atoi(cday); Int_t jn = atoi(cfno); Int_t inter = 0; char *zy; char *zm; char *zd; char *zn = ""; char *zyn; Int_t jyn = 0; Int_t goodthis = 0; Int_t onegood = 0; Long64_t ncalibs; TTree *tr; pamela::CalibCalPedEvent *ce = 0; pamela::PscuHeader *cph = 0; pamela::EventHeader *ceh = 0; if ( jn > 1 ) { jn--; } else { jn = 999; if ( jd > 1 ) { jd--; } else { if ( jm > 1 ) { jm--; } else { jm = 12; if ( jy > 3 ) { jy--; } else { return(1); }; }; if ( jm==4 || jm==6 || jm==9 || jm==11 ){ jd = 30; } else { jd = 31; }; if ( jm==2 ) jd = 29; }; }; Int_t gday = 0; Int_t retvar = 1; stringstream nyd; while ( jy>3 ) { while ( jm>0 ){ while ( jd > 0 ){ while ( jn > 0 ){ retvar--; nyear.str(""); nyear << jy; if ( jy<10 ){ zy = "0"; } else { zy = ""; }; nmonth.str(""); nmonth << jm; if ( jm<10 ){ zm = "0"; } else { zm = ""; }; nday.str(""); nday << jd; if ( jd<10 ){ zd = "0"; } else { zd = ""; }; nfno.str(""); nfno << jn; if ( jn>99 ) zn = ""; if ( jn<100 && jn>9 ) zn = "0"; if ( jn<10 ) zn = "00"; jyn = 0; goodthis = 0; onegood = 0; while ( jyn < 100 ){ begwhile: nyd.str(""); nyd << jyn; if ( jyn<10 ){ zyn = "0"; } else { zyn = ""; }; stringstream nfile; nfile.str(""); if ( upper ){ nfile << "DW_" << zy; } else { nfile << "dw_" << zy; }; nfile << nyear.str().c_str(); nfile << zm; nfile << nmonth.str().c_str(); nfile << zd; nfile << nday.str().c_str() << "_"; nfile << zn; nfile << nfno.str().c_str(); nfile << zyn; nfile << nyd.str().c_str(); // TString fgfile = ""; stringcopy(fgfile,ffile,0,posiz); stringappend(fgfile,nfile.str().c_str()); TFile *calheadFile = 0; calheadFile = getEmiFile(fgfile, "CalibCalPed", "Header"); TFile *calcalibFile = 0; calcalibFile = getEmiFile(fgfile, "CalibCalPed","Event"); // if ( !calheadFile || !calcalibFile ) { if ( calheadFile ) calheadFile->Close(); if ( calcalibFile ) calcalibFile->Close(); if ( jyn > 0 && onegood ) { goodthis = 1; jyn--; goto begwhile; }; if ( jyn == 0 ) { onegood = 0; jyn = 100; goto salta; }; if ( jyn > 0 && !onegood ) { onegood = 0; jyn = 100; goto salta; }; goto salta; } else { onegood = 1; if ( !goodthis ) { calheadFile->Close(); calcalibFile->Close(); goto salta; }; }; //Takes the tree of the header file tr = (TTree*)calheadFile->Get("Pscu"); tr->AddFriend("CalibCalPed", calcalibFile); tr->SetBranchAddress("Header", &ceh); tr->SetBranchAddress("Event", &ce); ncalibs = tr->GetEntries(); if ( ncalibs == 0 ) { jyn = 100; onegood = 0; goodthis = 0; calheadFile->Close(); calcalibFile->Close(); goto salta; }; inter = 0; for (Int_t c = 0; c < ncalibs; c++){ tr->GetEntry(c); cph = ceh->GetPscuHeader(); if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { inter++; }; }; if ( inter ){ jyn = 100; calheadFile->Close(); calcalibFile->Close(); return(retvar); }; calheadFile->Close(); calcalibFile->Close(); salta: jyn++; }; jn--; }; jn = 999; jd--; gday++; if ( gday>1 ) return(1); }; jm--; if ( jm==4 || jm==6 || jm==9 || jm==11 ){ jd = 30; } else { jd = 31; }; if ( jm==2 ) jd = 29; }; jm = 12; jy--; }; return(retvar); } TString whatnamewith(TString ffile, Int_t retvar){ Int_t debug = 0; const char *file = ffile; stringstream fil; fil.str(""); fil << file ; if ( debug ) printf("ffile %s retvar %i \n",file,retvar); Int_t posiz = fil.str().find("dw_"); Int_t upper = 0; if ( posiz == -1 ) { posiz = fil.str().find("DW_"); upper = 1; }; if ( posiz == -1 ) return("wrong"); // // printf("posiz %i \n",posiz); TString nomefile = getFilename(ffile); const char *ufnome = nomefile; stringstream nyear; stringstream nmonth; stringstream nday; stringstream nfno; TString year; TString month; TString day; TString fno; stringcopy(year, ufnome, 3, 5); stringcopy(month, ufnome, 5, 7); stringcopy(day, ufnome, 7, 9); stringcopy(fno, ufnome, 10, 13); const char *cyear = year; const char *cmonth = month; const char *cday = day; const char *cfno = fno; Int_t jy = atoi(cyear); Int_t jm = atoi(cmonth); Int_t jd = atoi(cday); Int_t jn = atoi(cfno); Int_t jyn = 0; Int_t goodthis = 0; Int_t onegood = 0; char *zy; char *zm; char *zd; char *zn = ""; char *zyn; Long64_t ncalibs; TString fgfile; TTree *tr; TFile *calheadFile = 0; TFile *calcalibFile = 0; if ( jn > 1 ) { jn--; } else { jn = 999; if ( jd > 1 ) { jd--; } else { if ( jm > 1 ) { jm--; } else { jm = 12; if ( jy > 3 ) { jy--; } else { return("wrong"); }; }; if ( jm==4 || jm==6 || jm==9 || jm==11 ){ jd = 30; } else { jd = 31; }; if ( jm==2 ) jd = 29; }; }; Int_t retval = 1; Int_t gday = 0; stringstream nyd; while ( jy>3 ) { while ( jm>0 ){ while ( jd > 0 ){ while ( jn > 0 ){ retval--; if ( retval < retvar ) return("wrong"); nyear.str(""); nyear << jy; if ( jy<10 ){ zy = "0"; } else { zy = ""; }; nmonth.str(""); nmonth << jm; if ( jm<10 ){ zm = "0"; } else { zm = ""; }; nday.str(""); nday << jd; if ( jd<10 ){ zd = "0"; } else { zd = ""; }; nfno.str(""); nfno << jn; if ( jn>99 ) zn = ""; if ( jn<100 && jn>9 ) zn = "0"; if ( jn<10 ) zn = "00"; jyn = 0; goodthis = 0; onegood = 0; while ( jyn < 100 ){ begwhile: nyd.str(""); nyd << jyn; if ( jyn<10 ){ zyn = "0"; } else { zyn = ""; }; stringstream nfile; nfile.str(""); if ( upper ){ nfile << "DW_" << zy; } else { nfile << "dw_" << zy; }; nfile << nyear.str().c_str(); nfile << zm; nfile << nmonth.str().c_str(); nfile << zd; nfile << nday.str().c_str() << "_"; nfile << zn; nfile << nfno.str().c_str(); nfile << zyn; nfile << nyd.str().c_str(); // fgfile = ""; stringcopy(fgfile,ffile,0,posiz); stringappend(fgfile,nfile.str().c_str()); const char *tosee = fgfile; if ( debug ) printf("file: %s retval %i retvar %i \n",tosee,retval,retvar); calheadFile = getEmiFile(fgfile, "CalibCalPed", "Header"); calcalibFile = getEmiFile(fgfile, "CalibCalPed","Event"); // if ( !calheadFile || !calcalibFile ) { if ( calheadFile ) calheadFile->Close(); if ( calcalibFile ) calcalibFile->Close(); if ( jyn > 0 && onegood ) { goodthis = 1; jyn--; goto begwhile; }; if ( jyn == 0 ) { onegood = 0; jyn = 100; goto salta; }; if ( jyn > 0 && !onegood ) { onegood = 0; jyn = 100; goto salta; }; goto salta; } else { onegood = 1; if ( !goodthis ){ calheadFile->Close(); calcalibFile->Close(); goto salta; }; }; tr = (TTree*)calheadFile->Get("Pscu"); tr->AddFriend("CalibCalPed", calcalibFile); ncalibs = tr->GetEntries(); if ( ncalibs == 0 ) { jyn = 100; onegood = 0; goodthis = 0; calheadFile->Close(); calcalibFile->Close(); goto salta; }; calheadFile->Close(); calcalibFile->Close(); if ( retval == retvar ) return(fgfile); salta: jyn++; }; jn--; }; jn = 999; jd--; gday++; if ( gday>7 ) return("wrong"); }; jm--; if ( jm==4 || jm==6 || jm==9 || jm==11 ){ jd = 30; } else { jd = 31; }; if ( jm==2 ) jd = 29; }; jm = 12; jy--; }; return("wrong"); } char *getLEVname(TString filename, TString detc, TString numb, TString frame="root"){ // char *file; stringstream file; file.str(""); const char *num = numb; const char *det = detc; const char *fra = frame; string fil = (const char *)filename; Int_t posiz = fil.find("dw_"); if ( posiz == -1 ) posiz = fil.find("DW_"); if ( posiz == -1 ) return(0); Int_t spos = posiz; Int_t epos = posiz+15; TString tmptempf; stringcopy(tmptempf,filename,spos,epos); const char *tempf = tmptempf; file << tempf << ".Physics.Level"; file << num << "."; file << det << ".Event."; file << fra; const char *rfile = file.str().c_str(); return((char*)rfile); }; char *getfileLEVname(TString filename, TString detc, TString numb, TString frame="root"){ // char *file; stringstream file; const char *num = numb; const char *det = detc; const char *fra = frame; string fil = (const char *)filename; Int_t posiz = fil.find("dw_"); if ( posiz == -1 ) posiz = fil.find("DW_"); if ( posiz == -1 ) return(0); Int_t spos = posiz; Int_t epos = posiz+13; TString tmptempf; stringcopy(tmptempf,filename,spos,epos); const char *tempf = tmptempf; file.str(""); file << tempf << ".Physics.Level"; // file << "00.Physics.Level"; file << num << "."; file << det << ".Event."; file << fra; const char *rfile = file.str().c_str(); return((char*)rfile); }; void OLDCaloFindCalibs(TString &filename, Calib & calib){ for (Int_t s = 0; s < 4; s++){ for (Int_t d = 1; d<50; d++){ calib.ttime[s][d] = 0; if ( d < 49 ) calib.time[s][d] = 0; }; }; TFile *calheadFile = 0; calheadFile = getEmiFile(filename, "CalibCalPed", "Header"); TFile *calcalibFile = 0; calcalibFile = getEmiFile(filename, "CalibCalPed","Event"); // if ( !calheadFile || !calcalibFile ) { if ( calheadFile ) calheadFile->Close(); if ( calcalibFile ) calcalibFile->Close(); printf("No calibration header file! Exiting...\n"); return; }; //Takes the tree of the header file TTree *tr = (TTree*)calheadFile->Get("Pscu"); tr->AddFriend("CalibCalPed", calcalibFile); pamela::CalibCalPedEvent *ce = 0; pamela::PscuHeader *cph = 0; pamela::EventHeader *ceh = 0; tr->SetBranchAddress("Header", &ceh); tr->SetBranchAddress("Event", &ce); Long64_t ncalibs = tr->GetEntries(); Int_t inter; for (Int_t s = 0; s < 4; s++){ for (Int_t d = 1; d<50; d++){ calib.ttime[s][d] = 0; if ( d < 49 ) calib.time[s][d] = 0; }; }; for (Int_t s = 0; s < 4; s++){ inter = 0; for (Int_t c = 0; c < ncalibs; c++){ tr->GetEntry(c); cph = ceh->GetPscuHeader(); calib.ttime[s][inter] = 0; if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { calib.ttime[s][inter] = cph->GetOrbitalTime(); inter++; } else { if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) { printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]); }; }; }; if ( inter == 0 ){ printf(" ERROR: no suitable calibration for section %i in this file!\n",s); }; for (Int_t d = 1; d<50; d++){ if ( calib.ttime[s][d] != 0 ) { calib.time[s][d-1] = calib.ttime[s][d-1] + (int)((calib.ttime[s][d] - calib.ttime[s][d-1])/2.); } else { if ( d == 1 ) { calib.time[s][d-1] = 0; }; }; }; }; calheadFile->Close(); calcalibFile->Close(); } void CaloFindCalibs(TString &filename, TString &calibfile, Int_t &wused, Calib & calib){ // // First of all check if the file has a monotone growing OBT time function. // Int_t debug = 0; Int_t obtjump = 0; Int_t evjump[50]; Int_t OBT = 0; Int_t OOBT = 0; Int_t OBT1 = 0; Int_t LOBT = 0; pamela::PscuHeader *ph = 0; pamela::EventHeader *eh = 0; TFile *headerF = 0; headerF = emigetFile(filename, "Physics", "Header"); if ( !headerF ) { printf(" No physics file! \n "); return; }; TTree *ctr = (TTree*)headerF->Get("Pscu"); Long64_t nevents = ctr->GetEntries(); if ( nevents == 0 ) { printf(" The file is empty! \n "); return; }; for (Int_t i = 0; i < nevents; i++){ ctr->SetBranchAddress("Header", &eh); ctr->GetEntry(i); ph = eh->GetPscuHeader(); OBT = ph->GetOrbitalTime(); if ( !i ) OBT1 = OBT; if ( i == nevents-1 ) LOBT = OBT; if ( OOBT > OBT ) { evjump[obtjump] = i; obtjump++; }; OOBT = OBT; }; printf("\n Obtjump = %i - FIRST OBT %i - LAST OBT %i \n\n",obtjump,OBT1,LOBT); // TTree *tr; pamela::CalibCalPedEvent *ce = 0; pamela::PscuHeader *cph = 0; pamela::EventHeader *ceh = 0; Long64_t ncalibs; // for (Int_t s = 0; s < 4; s++){ for (Int_t d = 1; d<50; d++){ calib.ttime[s][d] = 0; if ( d < 49 ) calib.time[s][d] = 0; }; }; for (Int_t s = 0; s < 4; s++){ // Int_t firstlook = 0; repeatsearch: // Int_t inter = 0; TFile *calheadFile = 0; calheadFile = getEmiFile(filename, "CalibCalPed", "Header"); TFile *calcalibFile = 0; calcalibFile = getEmiFile(filename, "CalibCalPed","Event"); // if ( !calheadFile || !calcalibFile ) { if ( calheadFile ) calheadFile->Close(); if ( calcalibFile ) calcalibFile->Close(); printf(" No calibration in this file! \n"); if ( !firstlook ){ wused = 2; filename = calibfile; firstlook = 1; goto repeatsearch; } else { wused = 0; goto jte; }; }; wused = 1; //Takes the tree of the header file tr = (TTree*)calheadFile->Get("Pscu"); tr->AddFriend("CalibCalPed", calcalibFile); tr->SetBranchAddress("Header", &ceh); tr->SetBranchAddress("Event", &ce); // ncalibs = tr->GetEntries(); calib.obtjump = 0; inter = 0; for (Int_t c = 0; c < ncalibs; c++){ tr->GetEntry(c); cph = ceh->GetPscuHeader(); calib.ttime[s][inter] = 0; if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { calib.ttime[s][inter] = cph->GetOrbitalTime(); inter++; } else { if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) { printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]); }; }; }; jte: if ( inter == 0 ){ printf(" WARNING: no suitable calibration for section %i in this file!\n",s); printf(" I WILL SEARCH IN PREVIOUS FILES\n"); }; if ( inter > 50 ){ printf(" WARNING: cannot handle more than 50 calibrations for file!\n"); printf(" I WILL SEARCH IN PREVIOUS FILES\n"); inter = 0; }; if ( obtjump ){ calib.obtjump = 1; printf(" ERROR: jump in OBT! CPU restarded? \n"); printf(" WARNING: I will use only the first calibration (if any)\n"); }; // Int_t lastcal = 0; Int_t of = 0; if ( OBT1 < calib.ttime[s][0] || !inter ) { // // here we must look for a calibration in the previous file... // TString ffile = ""; stringcopy(ffile,filename); Int_t fcodep = fetchpreviousfile(ffile,s); if ( fcodep == 1 ) { if ( inter ){ of = -3; printf(" WARNING: Problems to find a good calibration for section %i \n",s); } else { of = -2; printf(" ERROR: I was not able to find any good calibration for section %i \n",s); }; } else { stringcopy(ffile,filename); TString pfile = whatnamewith(ffile,fcodep); struct Calib tcal; OLDCaloFindCalibs(pfile,tcal); for ( Int_t d = 0; d<51; d++ ){ if ( tcal.ttime[s][d] == 0 ){ lastcal = d-1; break; }; }; if ( !inter ){ calib.ttime[s][0] = tcal.ttime[s][lastcal]; calib.time[s][0] = OBT1; calib.time[s][1] = LOBT; calib.fcode[s][0] = fcodep; if ( debug ) printf("boh! %i \n",fcodep); if ( debug ) getchar(); of = -1; } else { if ( obtjump ){ if ( debug ) printf("boh2! %i \n",fcodep); if ( debug ) getchar(); inter = 1; calib.time[s][0] = OBT1; calib.time[s][1] = LOBT; calib.fcode[s][0] = 10; of = -5; for (Int_t d = 1; d<50; d++){ calib.ttime[s][d] = 0; if ( d>1 && d<49 ) calib.time[s][d] = 0; }; } else { if ( lastcal >= 0 && tcal.ttime[s][lastcal]0 ){ for (Int_t i = inter+1; i>0; i--){ calib.ttime[s][i] = calib.ttime[s][i-1]; }; calib.time[s][0] = OBT1; calib.ttime[s][0] = tcal.ttime[s][lastcal]; calib.fcode[s][0] = fcodep; of = 2; } else { calib.time[s][0] = OBT1; calib.fcode[s][0] = 10; of = 1; }; }; }; }; } else { of = 0; }; if ( debug ) printf("boh3! %i \n",of); Int_t of1 = 0; if ( inter && of != 0 && of != 2 && of != -2 && of != -5) { for (Int_t i = inter+1; i>0; i--){ calib.ttime[s][i] = calib.ttime[s][i-1]; }; of1 = 0; }; if ( of == -2 ) of = -1; if ( of == -3 ) { of = 1; calib.fcode[s][0] = 10; }; // if ( of != -1 && of != -5 ){ if ( of == 2 ) of = 1; for (Int_t d = 0; dClose(); if ( calcalibFile ) calcalibFile->Close(); }; } void Calo1stcalib(TString filename, Calib & calib, Int_t b[4]){ Float_t estrip[2][22][96]; // // this is the value of the mip for each strip. To be changed when we will have the real values // for (Int_t s=0; s<4;s++){ for (Int_t d = 0; d<50; d++){ calib.ttime[s][d] = 0 ; if ( d < 49 ) calib.time[s][d] = 0 ; }; }; for (Int_t m = 0; m < 2 ; m++ ){ for (Int_t k = 0; k < 22; k++ ){ for (Int_t l = 0; l < 96; l++ ){ calib.calped[m][k][l] = 0. ; estrip[m][k][l] = 0.; }; }; } // // first of all find the calibrations in the file // OLDCaloFindCalibs(filename, calib); // // print on the screen the results: // printf(" ---------------------------------------------------------- \n \n"); Int_t calibex = 0; for (Int_t s=0; s<4;s++){ Int_t stop = 0; for (Int_t d = 0; d<48; d++){ if ( calib.ttime[s][d] != 0 ) calibex++; if ( calib.time[s][0] != 0 ){ if ( d == 0 ) printf(" Section %i from time 0 to time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d]); if ( calib.time[s][d+1] != 0 ) { printf(" Section %i from time %i to time %i use calibration at time %i \n",s,calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d+1]); } else { if ( !stop ){ printf(" Section %i from time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d+1]); stop = 1; }; }; } else { if ( calib.ttime[s][d] != 0 ) printf(" Section %i from time 0 use calibration at time %i \n",s,calib.ttime[s][d]); }; }; printf("\n"); }; printf(" ---------------------------------------------------------- \n"); if ( calibex < 1 ) { printf("No full calibration data in this file, sorry!\n"); } else { // // calibrate before starting // for (Int_t s = 0; s < 4; s++){ b[s]=0; CaloPede(filename,s,calib.ttime[s][b[s]],calib); }; }; } void ColorMIP(Float_t mip, int& colo){ // printf("mip = %f \n",mip); if ( colo > 0 ){ colo = 10; if ( mip > 0.7 ) colo = 38; if ( mip > 2. ) colo = 4; if ( mip > 10. ) colo = 3; if ( mip > 100. ) colo = 2; if ( mip > 500. ) colo = 6; } else { colo = 10; if ( mip > 0.7 ) colo = 17; if ( mip > 2. ) colo = 15; if ( mip > 10. ) colo = 14; if ( mip > 100. ) colo = 13; if ( mip > 500. ) colo = 12; }; } void ColorTOFMIP(Float_t mip, int& colo){ // printf("mip = %f \n",mip); if ( colo > 0 ){ colo = 10; if ( mip > 0. ) colo = 38; if ( mip > 2. ) colo = 4; if ( mip > 10. ) colo = 3; if ( mip > 100. ) colo = 2; if ( mip > 500. ) colo = 6; } else { colo = 10; if ( mip > 0. ) colo = 17; if ( mip > 2. ) colo = 15; if ( mip > 10. ) colo = 14; if ( mip > 100. ) colo = 13; if ( mip > 500. ) colo = 12; }; } Int_t WhatToDo(int& i, int& doflag, int& jumpto, Long64_t nevents, Float_t lastevno, Float_t firstevno, const char *file, TString outDir, TString figty, TCanvas &figure) { char *bw; if ( jumpto == -1 ){ bw = "bw"; } else { bw = ""; }; jumpto = 0; char tellme[256]; stringstream input; char tellme2[256]; stringstream input2; char tellme3[256]; stringstream input3; stringstream out; out.str("x"); input2.str("zzzzzzzzzzzzzz"); input3.str("z"); input << out.str().c_str(); stringstream stc; stringstream stc2; while ( !strcmp(input.str().c_str(),out.str().c_str()) ) { printf("\nPress to continue, b to go backward, j to jump to a\nparticular event number, p to save the figure, q to quit: \n"); cin.getline(tellme,256); // input.str(""); input << tellme; out << "y"; // stc.str("b"); // if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { if ( i > 0 ) { printf("WARNING: going backward!\n\n"); doflag = 2; } else { printf("This is the first event, you can't go backward! \n"); out.str(""); out << input.str().c_str(); }; }; stc.str("j"); if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { printf("\n Do you want to jump to a progressive number [p] or to an event number [e]? "); cin.getline(tellme3,256); input3.str(""); input3 << tellme3; stc.str("p"); if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) { printf("\n Enter the progressive number you want to jump to: "); cin.getline(tellme2,256); input2.str(""); input2 << tellme2; Int_t j; j = atoi(input2.str().c_str()); if ( j < 1 || j > nevents+1 ) { printf("\n You can choose between 1 and %i \n",(int)nevents+1); out.str(""); out << input.str().c_str(); } else { printf("\n Jumping to progressive number %i\n\n",j); i = j-2; }; }; stc.str("e"); if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) { printf("\n Enter the event number you want to jump to: "); cin.getline(tellme2,256); input2.str(""); input2 << tellme2; Int_t j; j = atoi(input2.str().c_str()); if ( j < firstevno || j > lastevno ) { printf("\n You can choose between %i and %i \n",(int)firstevno,(int)lastevno); out.str(""); out << input.str().c_str(); } else { printf("\n Jumping to event number %i\n\n",j); jumpto = j; i = -1; }; }; stc.str("e"); stc2.str("p"); if ( strcmp(input3.str().c_str(),stc.str().c_str()) && strcmp(input3.str().c_str(),stc2.str().c_str()) ) { printf(" You must type or \"p\" or \"e\"\n"); out.str(""); out << input.str().c_str(); }; }; // stc.str("q"); if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { printf("Exiting...\n"); return(1); }; // stc.str("p"); if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) { printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n"); cin.getline(tellme2,256); input2.str(""); input2 << tellme2; out.str(""); out << input.str().c_str(); // TString filename = file; const string fil = (const char*)filename; Int_t posiz = fil.find("dw_"); if ( posiz == -1 ) posiz = fil.find("DW_"); Int_t posiz2 = posiz+13; TString file2; stringcopy(file2,filename,posiz,posiz2); const char *figrec = file2; const char *outdir = outDir; stringstream figsave; const char *ty = figty; figsave.str(""); figsave << outdir << "/"; figsave << ty << "_"; figsave << (i+1) << "_"; figsave << figrec; figsave << bw << "."; figsave << input2.str().c_str(); figure.SaveAs(figsave.str().c_str()); printf("\n"); }; }; return(0); } void PrintFigure(TString filename, TString outDir, TString figty, TString saveas, TCanvas& figure) { const string fil = (const char*)filename; Int_t posiz = fil.find("dw_"); if ( posiz == -1 ) posiz = fil.find("DW_"); Int_t posiz2 = posiz+13; TString file2; stringcopy(file2,filename,posiz,posiz2); // const char *figrec = file2; stringstream figsave; figsave.str(""); figsave << outDir.Data() << "/"; figsave << figrec << "_"; figsave << figty.Data() << "."; figsave << saveas.Data(); figure.SaveAs(figsave.str().c_str()); return; } Double_t langaufun(Double_t *x, Double_t *par) { // Numeric constants Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) Double_t mpshift = -0.22278298; // Landau maximum location // Control constants Double_t np = 100.0; // number of convolution steps Double_t sc = 3.; // convolution extends to +-sc Gaussian sigmas // Variables Double_t xx; Double_t mpc; Double_t fland; Double_t sum = 0.0; Double_t xlow,xupp; Double_t step; Double_t i; // MP shift correction mpc = par[1] - mpshift * par[0]; // Range of convolution integral xlow = x[0] - sc * par[3]; xupp = x[0] + sc * par[3]; step = (xupp-xlow) / np; // Convolution integral of Landau and Gaussian by sum for(i=1.0; i<=np/2; i++) { xx = xlow + (i-.5) * step; fland = TMath::Landau(xx,mpc,par[0]) / par[0]; sum += fland * TMath::Gaus(x[0],xx,par[3]); xx = xupp - (i-.5) * step; fland = TMath::Landau(xx,mpc,par[0]) / par[0]; sum += fland * TMath::Gaus(x[0],xx,par[3]); }; return (par[2] * step * sum * invsq2pi / par[3]); } TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF, TString drw="QR0"){ // Once again, here are the Landau * Gaussian parameters: // par[0]=Width (scale) parameter of Landau density // par[1]=Most Probable (MP, location) parameter of Landau density // par[2]=Total area (integral -inf to inf, normalization constant) // par[3]=Width (sigma) of convoluted Gaussian function // // Variables for langaufit call: // his histogram to fit // fitrange[2] lo and hi boundaries of fit range // startvalues[4] reasonable start values for the fit // parlimitslo[4] lower parameter limits // parlimitshi[4] upper parameter limits // fitparams[4] returns the final fit parameters // fiterrors[4] returns the final fit errors // ChiSqr returns the chi square // NDF returns ndf Int_t i; Char_t FunName[100]; sprintf(FunName,"Fitfcn_%s",his->GetName()); // TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); if (ffitold) delete ffitold; // TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); ffit->SetParameters(startvalues); ffit->SetParNames("Width","MP","Area","GSigma"); // for (i=0; i<4; i++) { ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); }; his->Fit(FunName,drw); // fit within specified range, use ParLimits, do not plot // ffit->GetParameters(fitparams); // obtain fit parameters for (i=0; i<4; i++) { fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors } ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2 NDF[0] = ffit->GetNDF(); // obtain ndf // return (ffit); // return fit function } //Double_t langaupro(Double_t *params, Double_t &maxx) { Double_t langaupro(Double_t *params) { // Seaches for the location (x value) at the maximum of the // Landau-Gaussian convolute and its full width at half-maximum. // Double_t p,x; Double_t step; Double_t l,lold; Int_t i = 0; Int_t MAXCALLS = 5000; // Search for maximum p = params[1] - 0.1 * params[0]; step = 0.05 * params[0]; lold = -2.0; l = -1.0; // while ( (l != lold) && (i < MAXCALLS) ) { while ( fabs(l-lold)>10e-12 && (i < MAXCALLS) ) { i++; lold = l; x = p + step; l = langaufun(&x,params); if ( l < lold ) step = -step/10.; p += step; // printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step); // printf(" l %.16g \n",l); }; // printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step); //printf(" l %.16g \n",l); if (i == MAXCALLS) return (-1.); // return(x); } void delay(int selftrig, int &del){ int i; i = 0; del = 300; // while ((i < 16) && ((selftrig & (0x8000>>i)) == 0)) { del += 50; i++; }; if (i > 15) del = 1100; } Double_t fitraw(Double_t *x, Double_t *par){ // Double_t fitval = par[3] + par[0]*TMath::Exp(-par[2]*(x[0]-par[1])); // return fitval; } int checkmerging(Chmerge & chmerge){ for (Int_t i = chmerge.fromflag; i<1000; i++){ if ( chmerge.pscuOBT == 0 && chmerge.pscucount == 0 ) return(1); if ( chmerge.pscuOBT == chmerge.lastOBT[i] && chmerge.pscucount == chmerge.lastcount[i] ){ chmerge.fromflag = i; printf("\n WARNING! %i event number %i at OBT %i already processed, skipped!\n",i,chmerge.pscucount,chmerge.pscuOBT); return(1); }; }; return(0); } int fillmerging(Chmerge & chmerge){ chmerge.lastOBT[chmerge.scanflag] = chmerge.pscuOBT; chmerge.lastcount[chmerge.scanflag] = chmerge.pscucount; chmerge.scanflag++; if ( chmerge.scanflag > 999 ) return(0); chmerge.lastOBT[chmerge.scanflag] = 0; chmerge.lastcount[chmerge.scanflag] = 0; return(0); } //short int checkifpretyodaversion(TString filename){ // string fil = (const char *)filename; // Int_t posiz = fil.find("dw_"); // if ( posiz == -1 ) posiz = fil.find("DW_"); // if ( posiz == -1 ) return (-1); // TString year; // stringcopy(year,filename,posiz+3,posiz+5); // const char *ye = year; // Int_t y = atoi(ye); // TString month; // stringcopy(month,filename,posiz+5,posiz+7); // const char *mo = month; // Int_t m = atoi(mo); // TString day; // stringcopy(day,filename,posiz+7,posiz+9); // const char *da = day; // Int_t d = atoi(da); // TString fno; // stringcopy(fno,filename,posiz+10,posiz+13); // const char *fn = fno; // Int_t f = atoi(fn); // if ( y >= 5 && m >= 5 && d < 15 ) return(1); // if ( y == 5 && m == 5 && d == 15 && f < 7 ) return(1); // return(0); //}