#define NPMTs 24 int MakeTrendT0( char *infile, int run, char* ocdbStorage="raw://") { gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); gSystem->Load("libCORRFW"); gSystem->Load("libTender"); gSystem->Load("libPWGPP"); char *outfile = "trending.root"; if(!infile) return -1; if(!outfile) return -1; TFile *f = TFile::Open(infile,"read"); if (!f) { printf("File %s not available\n", infile); return -1; } // LOAD HISTOGRAMS FROM QAresults.root TObjArray *fTzeroObject = (TObjArray*) f->Get("T0_Performance/QAT0chists"); TH1F* fTzeroORAplusORC =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORAplusORC"))->Clone("A"); TH1F* fResolution =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fResolution"))->Clone("B"); TH1F* fTzeroORA =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORA"))->Clone("C"); TH1F* fTzeroORC =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORC"))->Clone("D"); TH2F *fTimeVSAmplitude[NPMTs];//counting PMTs from 0 TH1D *fAmplitude[NPMTs]; TH1D *fTime[NPMTs]; //-----> add new histogram here // sigma fit of resolution, mean T0A, T0C and T0C&A, mean amplitude for each PMT double resolutionSigma = -9999; //dummy init double tzeroOrAPlusOrC = -9999; //dummy init double tzeroOrA = -9999; //dummy init double tzeroOrC = -9999; //dummy init double meanAmplitude[NPMTs]; //dummy init double meanTime[NPMTs]; //dummy init double timeDelayOCDB[NPMTs]; //dummy init //................................................................................. for(int ipmt=1; ipmt<=NPMTs; ipmt++){ //loop over all PMTs fTimeVSAmplitude[ipmt-1] =(TH2F*) ((TH2F*) fTzeroObject->FindObject(Form("fTimeVSAmplitude%d",ipmt)))->Clone(Form("E%d",ipmt)); int nbinsX = fTimeVSAmplitude[ipmt-1]->GetNbinsX(); int nbinsY = fTimeVSAmplitude[ipmt-1]->GetNbinsY(); //Amplitude fAmplitude[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionX(Form("fAmplitude%d",ipmt), 1, nbinsY); meanAmplitude[ipmt-1] = -9999; //dummy init if(fAmplitude[ipmt-1]->GetEntries()>0){ meanAmplitude[ipmt-1] = fAmplitude[ipmt-1]->GetMean(); } //Time fTime[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionY(Form("fTime%d",ipmt), 1, nbinsX); meanTime[ipmt-1] = -9999; //dummy init if(fTime[ipmt-1]->GetEntries()>0){ if(fTime[ipmt-1]->GetEntries()>20){ meanTime[ipmt-1] = GetParameterGaus((TH1F*) fTime[ipmt-1], 1); // Mean Time }else if(fTime[ipmt-1]->GetEntries()>0){ meanTime[ipmt-1] = fTime[ipmt-1]->GetMean(); } } } if(fResolution->GetEntries()>20){ resolutionSigma = GetParameterGaus(fResolution, 2); //gaussian sigma }else if(fResolution->GetEntries()>0){ resolutionSigma = fResolution->GetRMS(); //gaussian sigma } if(fTzeroORAplusORC->GetEntries()>20){ tzeroOrAPlusOrC = GetParameterGaus(fTzeroORAplusORC, 1); //gaussian mean }else if(fTzeroORAplusORC->GetEntries()>0){ tzeroOrAPlusOrC = fTzeroORAplusORC->GetMean(); } if(fTzeroORA->GetEntries()>20){ tzeroOrA = GetParameterGaus(fTzeroORA, 1); //gaussian mean }else if(fTzeroORA->GetEntries()>0){ tzeroOrA = fTzeroORA->GetMean(); } if(fTzeroORC->GetEntries()>20){ tzeroOrC = GetParameterGaus(fTzeroORC, 1); //gaussian mean }else if(fTzeroORC->GetEntries()>0){ tzeroOrC = fTzeroORC->GetMean(); //gaussian mean } //-----> analyze the new histogram here and set mean/sigma f->Close(); //-------------------- READ OCDB TIME DELAYS --------------------------- // Arguments: AliCDBManager* man = AliCDBManager::Instance(); man->SetDefaultStorage(ocdbStorage); man->SetRun(run); AliCDBEntry *entry = AliCDBManager::Instance()->Get("T0/Calib/TimeDelay"); AliT0CalibTimeEq *clb = (AliT0CalibTimeEq*)entry->GetObject(); for (Int_t i=0; iGetCFDvalue(i,0); } //--------------- write walues to the output ------------ TTreeSRedirector* pcstream = NULL; pcstream = new TTreeSRedirector(outfile); if (!pcstream) return; TFile *x = pcstream->GetFile(); x->cd(); TObjString runType; Int_t startTimeGRP=0; Int_t stopTimeGRP=0; Int_t time=0; Int_t duration=0; time = (startTimeGRP+stopTimeGRP)/2; duration = (stopTimeGRP-startTimeGRP); (*pcstream)<<"t0QA"<< "run="<Close(); delete pcstream; return; } //_____________________________________________________________________________ double GetParameterGaus(TH1F *histo, int whichParameter){ int maxBin = histo->GetMaximumBin(); double max = (histo->GetBinContent(maxBin-1) + histo->GetBinContent(maxBin) + histo->GetBinContent(maxBin+1))/3; double mean = histo->GetBinCenter(maxBin); //mean double lowfwhm = histo->GetBinCenter(histo->FindFirstBinAbove(max/2)); double highfwhm = histo->GetBinCenter(histo->FindLastBinAbove(max/2)); double sigma = (highfwhm - lowfwhm)/2.35482; //estimate fwhm FWHM = 2.35482*sigma TF1 *gaussfit = new TF1("gaussfit","gaus", mean - 4*sigma, mean + 4*sigma); // fit in +- 4 sigma window gaussfit->SetParameters(max, mean, sigma); if(whichParameter==2) histo->Fit(gaussfit,"RQNI"); else histo->Fit(gaussfit,"RQN"); double parValue = gaussfit->GetParameter(whichParameter); delete gaussfit; return parValue; }