[u/mrichter/AliRoot.git] / PMD / AliPMDCalibGain.cxx

//////////////////////////////////////////////////////////////////////////////
//
// * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
// *                                                                        *
// * Author: The ALICE Off-line Project.                                    *
// * Contributors are mentioned in the code where appropriate.  	    *
// *									    *
// * Permission to use, copy, modify and distribute this software and its   *
// * documentation strictly for non-commercial purposes is hereby granted   *
// * without fee, provided that the above copyright notice appears in all   *
// * copies and that both the copyright notice and this permission notice   *
// * appear in the supporting documentation. The authors make no claims     *
// * about the suitability of this software for any purpose. It is	    *
// * provided "as is" without express or implied warranty.		    *
// **************************************************************************/
//
//////////////////////////////////////////////////////////////////////////////

#include "Riostream.h"
#include "TF1.h"
#include "TFile.h"
#include "TObjString.h"
#include "TROOT.h"
#include "TClonesArray.h"
#include "TH1F.h"
#include "TObjArray.h"
#include "TTree.h"
#include "TMath.h"

// --- Standard library ---

// --- AliRoot header files ---
#include "AliDAQ.h"
#include "AliLog.h"
#include "AliRawReader.h"
#include "AliPMDRawStream.h"
#include "AliPMDddldata.h"
#include "AliPMDCalibGain.h"

ClassImp(AliPMDCalibGain)

AliPMDCalibGain::AliPMDCalibGain():
  TObject(),
  fpw(NULL)
{
  // Standard Constructor
  for(Int_t idet = 0; idet < kDet; idet++)
    {
      for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	{
	  fSMIso[idet][ismn]   = 0.;
	  fSMCount[idet][ismn] = 0.;
	  fCountSm[idet][ismn]=0.;
	  fTempnhit[idet][ismn]=0.;
	  fTempnhitSq[idet][ismn]=0.;
	  for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
	    {
	      for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
		{
		  fCellIso[idet][ismn][jrow][kcol]    = 0.;
		  fCellCount[idet][ismn][jrow][kcol]  = 0.;
		  fNhitCell[idet][ismn][jrow][kcol]   = 0.;
		  fPedMeanRMS[idet][ismn][jrow][kcol] = 0.;
		  fHotFlag[idet][ismn][jrow][kcol]    = 0.;
		  
		}
	    }
	}
    }
  
}
// ------------------------------------------------------------------------ //
AliPMDCalibGain::AliPMDCalibGain(const AliPMDCalibGain &pmdcalibgain):
  TObject(pmdcalibgain),
  fpw(NULL)
{
  for(Int_t idet = 0; idet < kDet; idet++)
    {
      for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	{
	  fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn] ;
	  fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn] ;
	  fCountSm[idet][ismn] = pmdcalibgain.fCountSm[idet][ismn];
	  fTempnhit[idet][ismn] = pmdcalibgain.fTempnhit[idet][ismn];
	  fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
	  for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
	    {
	      for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
		{
		  fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
		  fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
		  fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol];
		  fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
		  fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
		  
		}
	    }
	}
    }
  
}
// ------------------------------------------------------------------------ //
AliPMDCalibGain &AliPMDCalibGain::operator=(const AliPMDCalibGain &pmdcalibgain)
{
  if(this != &pmdcalibgain)
    {
      this->fpw = pmdcalibgain.fpw;
      for(Int_t idet = 0; idet < kDet; idet++)
	{
	  for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	    {
	      fSMIso[idet][ismn]      = pmdcalibgain.fSMIso[idet][ismn];
	      fSMCount[idet][ismn]    = pmdcalibgain.fSMCount[idet][ismn];
	      fCountSm[idet][ismn]    = pmdcalibgain.fCountSm[idet][ismn];
	      fTempnhit[idet][ismn]   = pmdcalibgain.fTempnhit[idet][ismn];
	      fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
	      for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
		{
		  for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
		    {
		      fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
		      fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
		      fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol]; //za 
		      fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
		      fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
		      
		    }
		}
	    }
	}
    }
  return *this;
}
// ------------------------------------------------------------------------ //
AliPMDCalibGain::~AliPMDCalibGain()
{
  // dtor
}

// ------------------------------------------------------------------------ //

Int_t AliPMDCalibGain::ExtractPedestal(const Char_t *rootFile)
{
  // Pedestal extraction from the PMD_PED.root file
  // To be called once at the beginning
  
  Int_t   det=0, sm=0, row=0, col=0;
  Float_t mean=0., rms=0.;
  
  TFile *pedfile = new TFile(rootFile);
  
  if(!pedfile)
    {
      printf("ERROR --- NO PEDESTAL (PMD_PED1.root) FILE IS FOUND --- STOP GAIN DA\n");
      return -3;
    }


  TTree *ped =(TTree*)pedfile->Get("ped");

  ped->SetBranchAddress("det",&det);
  ped->SetBranchAddress("sm",&sm);
  ped->SetBranchAddress("row",&row);
  ped->SetBranchAddress("col",&col);
  ped->SetBranchAddress("mean",&mean);
  ped->SetBranchAddress("rms",&rms);

  Int_t nentries = (Int_t)ped->GetEntries();

  for (Int_t ient = 0; ient < nentries; ient++)
    {
      ped->GetEntry(ient);
      fPedMeanRMS[det][sm][row][col] = mean + 3.*rms;
      //printf("Mean= %f, RMS= %f, PedMeanRMS=%f\n",mean,rms,fPedMeanRMS[det][sm][row][col]);

    }

  pedfile->Close();
  delete pedfile;
  pedfile = 0x0;

  return 1;
}
//------------------------------------------------------------------------------------------------

Int_t AliPMDCalibGain::ExtractHotChannel(const Char_t *rootFile)
{
  // HotChannel extraction from the PMD_HOT.root file
  // To be called once at the beginning

  Int_t   det=0, sm=0, row=0, col=0;
  Float_t flag=0.;

  TFile *hotmapfile = new TFile(rootFile);

  if(!hotmapfile)
    {
      printf(" NO HOTCHANNEL MAP (PMD_HOT.root) FILE IS FOUND \n");

      for (Int_t idet = 0; idet < kDet; idet++)
	{
	  for (Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	    {
	      for (Int_t irow = 0; irow < kMaxRow; irow++)
		{
		  for (Int_t icol = 0; icol < kMaxCol; icol++)
		    {
		      fHotFlag[idet][ismn][irow][icol] = 0.;
		    }
		}
	    }
	}
    }


  TTree *hot =(TTree*)hotmapfile->Get("hot");
  
  hot->SetBranchAddress("det",&det);
  hot->SetBranchAddress("sm",&sm);
  hot->SetBranchAddress("row",&row);
  hot->SetBranchAddress("col",&col);
  hot->SetBranchAddress("flag",&flag);
  
  Int_t nentries = (Int_t)hot->GetEntries();
  
  for (Int_t ient = 0; ient < nentries; ient++)
    {
      hot->GetEntry(ient);
      fHotFlag[det][sm][row][col] = flag;

      //printf(" HotFlag=%f\n",fHotFlag[det][sm][row][col]);
    }

  hotmapfile->Close();
  delete hotmapfile;
  hotmapfile = 0x0;
  
  return 1;
}


// ------------------------------------------------------------------------ //

void AliPMDCalibGain::ReadTempFile(const Char_t *tempFile)
{
  // Read the variables from the file
  
  ifstream intmpfile;
  intmpfile.open(tempFile);

  Int_t iddet = 0, issm = 0, irrow = 0, iccol = 0;
  Float_t smcount = 0., smiso = 0.;
  Float_t cellcount = 0., celliso = 0.;


  for (Int_t idet = 0; idet < kDet; idet++)
    {
      for (Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  intmpfile >> iddet >> issm >> smcount >> smiso;
	  fSMCount[idet][ism] = smcount;
	  fSMIso[idet][ism]   = smiso;
	}
    }

  for (Int_t idet = 0; idet < kDet; idet++)
    {
      for (Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  for (Int_t irow = 0; irow < kMaxRow; irow++)
	    {
	      for (Int_t icol = 0; icol < kMaxCol; icol++)
		{
		  intmpfile >> iddet >> issm >> irrow >> iccol
			    >> cellcount >> celliso;
		  fCellCount[idet][ism][irow][icol] = cellcount;
		  fCellIso[idet][ism][irow][icol]   = celliso;
		}
	    }
	}
    }
  

  intmpfile.close();

}
// ------------------------------------------------------------------------ //
void AliPMDCalibGain::WriteTempFile(const Char_t *tempFile)
{
  // Write the Temporary file if the required statics is not achieved


  /*
    Following variables to be written to a file
    fDetIso[idet] ;
    fSMIso[idet][ismn]; 
    fCellIso[idet][ismn][irow][icol]; 
    
    fDetCount[idet];
    fSMCount[idet][ismn];
    fCellCount[idet][ismn][irow][icol];
  */				  


  fpw = fopen(tempFile,"w+");
  
  for (Int_t idet = 0; idet < kDet; idet++)
    {
     // fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
    }
  
  for (Int_t idet = 0; idet < kDet; idet++)
    {
      for (Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
	}
    }
  
  for (Int_t idet = 0; idet < kDet; idet++)
    {
      for (Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  for (Int_t irow = 0; irow < kMaxRow; irow++)
	    {
	      for (Int_t icol = 0; icol < kMaxCol; icol++)
		{
		  fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
			  fCellCount[idet][ism][irow][icol],
			  fCellIso[idet][ism][irow][icol]);
		}
	    }
	}
    }
  
  fclose(fpw);
  
}

// ------------------------------------------------------------------------ //

Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader *rawReader, TObjArray *pmdddlcont)
{
  // Calculates the ADC of isolated cell
  
  const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
  const Int_t kCellNeighbour = 6;
  Int_t neibx[6] = {1,0,-1,-1,0,1};
  Int_t neiby[6] = {0,1,1,0,-1,-1};
  
  Int_t id1 = 0,jd1 = 0;  //neighbour row/col
  Int_t isocount = 0;     //number of neighbours with 0 signal

  Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
  
  for(Int_t idet = 0; idet < kDet; idet++)
    {
      for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
        {
	  for(Int_t irow = 0; irow < kMaxRow; irow++)
            {
              for(Int_t icol = 0; icol < kMaxCol; icol++)
                {
		  d1[idet][ismn][irow][icol] = 0.;
                }
	    }
	}
    }
  
  AliPMDRawStream rawStream(rawReader);

  Int_t iddl = -1;
  
  Int_t numberofDDLs = 0;

  while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
    numberofDDLs++;
    
    Int_t ientries = pmdddlcont->GetEntries();
    
    for (Int_t ient = 0; ient < ientries; ient++)
      {
	AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont->UncheckedAt(ient);
	
	Int_t idet = pmdddl->GetDetector();
	Int_t ismn = pmdddl->GetSMN();
	Int_t mcm = pmdddl->GetMCM();
	//Int_t ichno = pmdddl->GetChannel();
	Int_t irow = pmdddl->GetRow();
	Int_t icol = pmdddl->GetColumn();
	Int_t isig = pmdddl->GetSignal();
	
	// This is the protection not to crash the code 
	
	if(mcm == 0) continue;
	if (irow < 0 || icol < 0 || irow > 47 || icol > 95) continue;
	
	// Pedestal subtraction
	if(fHotFlag[idet][ismn][irow][icol] == 1.0) isig = 0;
	
	if (isig>0)
	  {
	    d1[idet][ismn][irow][icol] =
	      (Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
	    
	    //printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
	    
	    fNhitCell[idet][ismn][irow][icol]++;     //cell hit frequency
	  }
	

      }
    pmdddlcont->Delete();
    }
  
  for(Int_t idet=0; idet < kDet; idet++)
    {
      for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	{
	  for(Int_t irow = 0; irow < kMaxRow; irow++)
	    {
	      for(Int_t icol = 0; icol < kMaxCol; icol++)
		{
		  if(d1[idet][ismn][irow][icol] > 0)
		    {
		      isocount = 0;
		      for(Int_t ii = 0; ii < kCellNeighbour; ii++)
			{
			  id1 = irow + neibx[ii];
			  jd1 = icol + neiby[ii];
			  if (id1 < 0) id1 = 0;
			  if (id1 > kMaxRow-1) id1 = kMaxRow - 1;
			  if (jd1 < 0) jd1 = 0;
			  if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1;
			  if(d1[idet][ismn][id1][jd1] == 0)
			    {
			      isocount++;
			      if(isocount == kCellNeighbour)
				{
				  //fDetIso[idet] += d1[idet][ismn][irow][icol];
				  fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
				  fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
				  //fDetCount[idet]++;
				  fSMCount[idet][ismn]++;
				  fCellCount[idet][ismn][irow][icol]++;
				  
				}
			    }
			}  // neigh cell cond.
		    }     // d>0 cond.
		}
	    }
	}
    }
  
  for(Int_t idet=0; idet < kDet; idet++)
    {
      for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	{
	  for(Int_t irow = 0; irow < kMaxRow; irow++)
	    {
	      for(Int_t icol = 0; icol < kMaxCol; icol++)
		{
		  if(fNhitCell[idet][ismn][irow][icol]>0)
		    {
		      fCountSm[idet][ismn]  += 1;
		      fTempnhit[idet][ismn] += fNhitCell[idet][ismn][irow][icol];
		      fTempnhitSq[idet][ismn] += fNhitCell[idet][ismn][irow][icol]
			*fNhitCell[idet][ismn][irow][icol];
		    }
		}
	    }
	}
    }
  
  if (numberofDDLs < kDDL) return kFALSE;
  return kTRUE;
  
}
// ------------------------------------------------------------------------ //
void AliPMDCalibGain::Analyse(TTree *gaintree, TTree *meantree)
{
  // Calculates the mean
  Int_t   det = 0, sm = 0, row = 0, col = 0;
  Float_t gain = 0.;
  Float_t cellmean = 0.;

  Float_t modmean[2][24];

  for (Int_t idet=0; idet < 2; idet++)
    {
      for (Int_t ism = 0; ism < 24; ism++)
	{
	  modmean[idet][ism] = 0.;
	}
    }

  gaintree->Branch("det",&det,"det/I");
  gaintree->Branch("sm",&sm,"sm/I");
  gaintree->Branch("row",&row,"row/I");
  gaintree->Branch("col",&col,"col/I");
  gaintree->Branch("gain",&gain,"gain/F");
  
  for(Int_t idet = 0; idet < kDet; idet++)
    {
      for(Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  if (fSMCount[idet][ism] > 0)
	    modmean[idet][ism] = fSMIso[idet][ism]/fSMCount[idet][ism];
	  for(Int_t irow = 0; irow < kMaxRow; irow++)
	    {
	      for(Int_t icol = 0; icol < kMaxCol; icol++)
		{
		  if (fCellCount[idet][ism][irow][icol] > 0.)
		    {
		      cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
		    }
		  det      = idet;
		  sm       = ism;
		  row      = irow;
		  col      = icol;
		  if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
		    {
		      gain = cellmean/modmean[idet][ism];
		    }
		  else
		    {
		      gain = 0.;
		    }
		  //if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
		  gaintree->Fill();
		}
	    }
	}
    }

  Float_t smmean;

  // Writing each module mean value
  meantree->Branch("det",&det,"det/I");
  meantree->Branch("sm",&sm,"sm/I");
  meantree->Branch("smmean",&smmean,"row/F");
  
  for(Int_t idet = 0; idet < kDet; idet++)
    {
      for (Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  det    = idet;
	  sm     = ism;
	  smmean = modmean[idet][ism];
	  meantree->Fill();
	}
    }

}
// ------------------------------------------------------------------------ //
void AliPMDCalibGain::FindHotCell(TTree *hottree, Float_t xvar)
{
  // Calculates the mean
  Int_t   det = 0, sm = 0, row = 0, col = 0;
  Float_t flag = 0.;
  Float_t meannhit = 0.;
  Float_t meanSqnhit = 0.;
  Float_t sigmanhit = 0.,nhitcut = 0.;

  //Float_t xvar = 5.;

  hottree->Branch("det",&det,"det/I");
  hottree->Branch("sm",&sm,"sm/I");
  hottree->Branch("row",&row,"row/I");
  hottree->Branch("col",&col,"col/I");
  hottree->Branch("flag",&flag,"flag/F");
  
  for(Int_t idet = 0; idet < kDet; idet++)
    {
      for(Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  if (fCountSm[idet][ism]> 0)
	    {
	      meannhit   = fTempnhit[idet][ism]/fCountSm[idet][ism];
	      meanSqnhit = fTempnhitSq[idet][ism]/fCountSm[idet][ism];
	      sigmanhit  = sqrt(meanSqnhit-(meannhit*meannhit));
	      nhitcut    = meannhit + xvar*sigmanhit;

	      for(Int_t irow = 0; irow < kMaxRow; irow++)
		{
		  for(Int_t icol = 0; icol < kMaxCol; icol++)
		    {
		      det      = idet;
		      sm       = ism;
		      row      = irow;
		      col      = icol;
		      
		      if(fNhitCell[idet][ism][irow][icol] > nhitcut)
			{
			  flag  = 1.0;
			}
		      else
			{
			  flag = 0.;
			}
		      hottree->Fill();
		    }
		  
		}
	    }
	}
    }
}
Commit	Line	Data
4fc06b3a	1	//////////////////////////////////////////////////////////////////////////////
	2	//
	3	// * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	// * *
	5	// * Author: The ALICE Off-line Project. *
	6	// * Contributors are mentioned in the code where appropriate. *
	7	// * *
	8	// * Permission to use, copy, modify and distribute this software and its *
	9	// * documentation strictly for non-commercial purposes is hereby granted *
	10	// * without fee, provided that the above copyright notice appears in all *
	11	// * copies and that both the copyright notice and this permission notice *
	12	// * appear in the supporting documentation. The authors make no claims *
	13	// * about the suitability of this software for any purpose. It is *
	14	// * provided "as is" without express or implied warranty. *
	15	// **************************************************************************/
	16	//
	17	//////////////////////////////////////////////////////////////////////////////
	18
4e2efa25	19	#include "Riostream.h"
4fc06b3a	20	#include "TF1.h"
	21	#include "TFile.h"
	22	#include "TObjString.h"
	23	#include "TROOT.h"
	24	#include "TClonesArray.h"
	25	#include "TH1F.h"
	26	#include "TObjArray.h"
	27	#include "TTree.h"
7c11ad1f	28	#include "TMath.h"
4fc06b3a	29
	30	// --- Standard library ---
	31
	32	// --- AliRoot header files ---
	33	#include "AliDAQ.h"
	34	#include "AliLog.h"
	35	#include "AliRawReader.h"
	36	#include "AliPMDRawStream.h"
	37	#include "AliPMDddldata.h"
	38	#include "AliPMDCalibGain.h"
	39
	40	ClassImp(AliPMDCalibGain)
	41
1da50520	42	AliPMDCalibGain::AliPMDCalibGain():
	43	TObject(),
	44	fpw(NULL)
4fc06b3a	45	{
4fc06b3a	46	// Standard Constructor
7c11ad1f	47	for(Int_t idet = 0; idet < kDet; idet++)
4fc06b3a	48	{
7c11ad1f	49	for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
4fc06b3a	50	{
7c11ad1f	51	fSMIso[idet][ismn] = 0.;
	52	fSMCount[idet][ismn] = 0.;
	53	fCountSm[idet][ismn]=0.;
	54	fTempnhit[idet][ismn]=0.;
	55	fTempnhitSq[idet][ismn]=0.;
	56	for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
4fc06b3a	57	{
7c11ad1f	58	for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
4fc06b3a	59	{
7c11ad1f	60	fCellIso[idet][ismn][jrow][kcol] = 0.;
	61	fCellCount[idet][ismn][jrow][kcol] = 0.;
	62	fNhitCell[idet][ismn][jrow][kcol] = 0.;
	63	fPedMeanRMS[idet][ismn][jrow][kcol] = 0.;
	64	fHotFlag[idet][ismn][jrow][kcol] = 0.;
	65
4fc06b3a	66	}
	67	}
	68	}
	69	}
7c11ad1f	70
4fc06b3a	71	}
	72	// ------------------------------------------------------------------------ //
	73	AliPMDCalibGain::AliPMDCalibGain(const AliPMDCalibGain &pmdcalibgain):
1da50520	74	TObject(pmdcalibgain),
1da50520	75	fpw(NULL)
4fc06b3a	76	{
7c11ad1f	77	for(Int_t idet = 0; idet < kDet; idet++)
4fc06b3a	78	{
7c11ad1f	79	for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
4fc06b3a	80	{
7c11ad1f	81	fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn] ;
	82	fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn] ;
	83	fCountSm[idet][ismn] = pmdcalibgain.fCountSm[idet][ismn];
	84	fTempnhit[idet][ismn] = pmdcalibgain.fTempnhit[idet][ismn];
	85	fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
	86	for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
4fc06b3a	87	{
7c11ad1f	88	for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
4fc06b3a	89	{
7c11ad1f	90	fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
	91	fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
	92	fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol];
	93	fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
	94	fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
	95
4fc06b3a	96	}
	97	}
	98	}
	99	}
7c11ad1f	100
4fc06b3a	101	}
	102	// ------------------------------------------------------------------------ //
	103	AliPMDCalibGain &AliPMDCalibGain::operator=(const AliPMDCalibGain &pmdcalibgain)
	104	{
7c11ad1f	105	if(this != &pmdcalibgain)
4fc06b3a	106	{
1da50520	107	this->fpw = pmdcalibgain.fpw;
7c11ad1f	108	for(Int_t idet = 0; idet < kDet; idet++)
4fc06b3a	109	{
7c11ad1f	110	for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
4fc06b3a	111	{
7c11ad1f	112	fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn];
	113	fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn];
	114	fCountSm[idet][ismn] = pmdcalibgain.fCountSm[idet][ismn];
	115	fTempnhit[idet][ismn] = pmdcalibgain.fTempnhit[idet][ismn];
	116	fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
	117	for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
4fc06b3a	118	{
7c11ad1f	119	for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
	120	{
	121	fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
	122	fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
	123	fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol]; //za
	124	fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
	125	fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
	126
4fc06b3a	127	}
	128	}
	129	}
	130	}
	131	}
7c11ad1f	132	return *this;
4fc06b3a	133	}
	134	// ------------------------------------------------------------------------ //
	135	AliPMDCalibGain::~AliPMDCalibGain()
	136	{
7c11ad1f	137	// dtor
1da50520	138	}
	139
	140	// ------------------------------------------------------------------------ //
	141
d31b0352	142	Int_t AliPMDCalibGain::ExtractPedestal(const Char_t *rootFile)
1da50520	143	{
	144	// Pedestal extraction from the PMD_PED.root file
	145	// To be called once at the beginning
7c11ad1f	146
be8b7039	147	Int_t det=0, sm=0, row=0, col=0;
be8b7039	148	Float_t mean=0., rms=0.;
7c11ad1f	149
d31b0352	150	TFile *pedfile = new TFile(rootFile);
7c11ad1f	151
1da50520	152	if(!pedfile)
	153	{
	154	printf("ERROR --- NO PEDESTAL (PMD_PED1.root) FILE IS FOUND --- STOP GAIN DA\n");
	155	return -3;
	156	}
	157
	158
	159	TTree ped =(TTree)pedfile->Get("ped");
	160
	161	ped->SetBranchAddress("det",&det);
	162	ped->SetBranchAddress("sm",&sm);
	163	ped->SetBranchAddress("row",&row);
	164	ped->SetBranchAddress("col",&col);
	165	ped->SetBranchAddress("mean",&mean);
	166	ped->SetBranchAddress("rms",&rms);
	167
	168	Int_t nentries = (Int_t)ped->GetEntries();
	169
	170	for (Int_t ient = 0; ient < nentries; ient++)
	171	{
	172	ped->GetEntry(ient);
	173	fPedMeanRMS[det][sm][row][col] = mean + 3.*rms;
	174	//printf("Mean= %f, RMS= %f, PedMeanRMS=%f\n",mean,rms,fPedMeanRMS[det][sm][row][col]);
	175
	176	}
	177
	178	pedfile->Close();
	179	delete pedfile;
	180	pedfile = 0x0;
	181
	182	return 1;
	183	}
7c11ad1f	184	//------------------------------------------------------------------------------------------------
	185
	186	Int_t AliPMDCalibGain::ExtractHotChannel(const Char_t *rootFile)
	187	{
	188	// HotChannel extraction from the PMD_HOT.root file
	189	// To be called once at the beginning
	190
be8b7039	191	Int_t det=0, sm=0, row=0, col=0;
be8b7039	192	Float_t flag=0.;
7c11ad1f	193
	194	TFile *hotmapfile = new TFile(rootFile);
	195
	196	if(!hotmapfile)
	197	{
637f01e6	198	printf(" NO HOTCHANNEL MAP (PMD_HOT.root) FILE IS FOUND \n");
	199
	200	for (Int_t idet = 0; idet < kDet; idet++)
	201	{
4e2efa25	202	for (Int_t ismn = 0; ismn < kMaxSMN; ismn++)
637f01e6	203	{
4e2efa25	204	for (Int_t irow = 0; irow < kMaxRow; irow++)
637f01e6	205	{
4e2efa25	206	for (Int_t icol = 0; icol < kMaxCol; icol++)
637f01e6	207	{
a7307087	208	fHotFlag[idet][ismn][irow][icol] = 0.;
637f01e6	209	}
	210	}
	211	}
	212	}
7c11ad1f	213	}
	214
	215
	216	TTree hot =(TTree)hotmapfile->Get("hot");
	217
	218	hot->SetBranchAddress("det",&det);
	219	hot->SetBranchAddress("sm",&sm);
	220	hot->SetBranchAddress("row",&row);
	221	hot->SetBranchAddress("col",&col);
	222	hot->SetBranchAddress("flag",&flag);
	223
	224	Int_t nentries = (Int_t)hot->GetEntries();
	225
	226	for (Int_t ient = 0; ient < nentries; ient++)
	227	{
	228	hot->GetEntry(ient);
	229	fHotFlag[det][sm][row][col] = flag;
	230
	231	//printf(" HotFlag=%f\n",fHotFlag[det][sm][row][col]);
	232	}
	233
	234	hotmapfile->Close();
	235	delete hotmapfile;
	236	hotmapfile = 0x0;
	237
	238	return 1;
	239	}
	240
	241
1da50520	242	// ------------------------------------------------------------------------ //
1da50520	243
d31b0352	244	void AliPMDCalibGain::ReadTempFile(const Char_t *tempFile)
1da50520	245	{
	246	// Read the variables from the file
	247
76942fe9	248	ifstream intmpfile;
	249	intmpfile.open(tempFile);
	250
44372324	251	Int_t iddet = 0, issm = 0, irrow = 0, iccol = 0;
be8b7039	252	Float_t smcount = 0., smiso = 0.;
be8b7039	253	Float_t cellcount = 0., celliso = 0.;
1da50520	254
1da50520	255
	256	for (Int_t idet = 0; idet < kDet; idet++)
	257	{
	258	for (Int_t ism = 0; ism < kMaxSMN; ism++)
	259	{
76942fe9	260	intmpfile >> iddet >> issm >> smcount >> smiso;
1da50520	261	fSMCount[idet][ism] = smcount;
	262	fSMIso[idet][ism] = smiso;
	263	}
	264	}
	265
	266	for (Int_t idet = 0; idet < kDet; idet++)
	267	{
	268	for (Int_t ism = 0; ism < kMaxSMN; ism++)
	269	{
	270	for (Int_t irow = 0; irow < kMaxRow; irow++)
	271	{
	272	for (Int_t icol = 0; icol < kMaxCol; icol++)
	273	{
76942fe9	274	intmpfile >> iddet >> issm >> irrow >> iccol
76942fe9	275	>> cellcount >> celliso;
1da50520	276	fCellCount[idet][ism][irow][icol] = cellcount;
	277	fCellIso[idet][ism][irow][icol] = celliso;
	278	}
	279	}
	280	}
	281	}
7c11ad1f	282
76942fe9	283
76942fe9	284	intmpfile.close();
1da50520	285
4fc06b3a	286	}
4fc06b3a	287	// ------------------------------------------------------------------------ //
d31b0352	288	void AliPMDCalibGain::WriteTempFile(const Char_t *tempFile)
1da50520	289	{
d31b0352	290	// Write the Temporary file if the required statics is not achieved
1da50520	291
1da50520	292
1da50520	293	/*
d31b0352	294	Following variables to be written to a file
1da50520	295	fDetIso[idet] ;
	296	fSMIso[idet][ismn];
	297	fCellIso[idet][ismn][irow][icol];
	298
	299	fDetCount[idet];
	300	fSMCount[idet][ismn];
	301	fCellCount[idet][ismn][irow][icol];
	302	*/
	303
d31b0352	304
d31b0352	305	fpw = fopen(tempFile,"w+");
7c11ad1f	306
1da50520	307	for (Int_t idet = 0; idet < kDet; idet++)
1da50520	308	{
7c11ad1f	309	// fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
1da50520	310	}
7c11ad1f	311
1da50520	312	for (Int_t idet = 0; idet < kDet; idet++)
	313	{
	314	for (Int_t ism = 0; ism < kMaxSMN; ism++)
	315	{
	316	fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
	317	}
	318	}
7c11ad1f	319
1da50520	320	for (Int_t idet = 0; idet < kDet; idet++)
	321	{
	322	for (Int_t ism = 0; ism < kMaxSMN; ism++)
	323	{
	324	for (Int_t irow = 0; irow < kMaxRow; irow++)
	325	{
	326	for (Int_t icol = 0; icol < kMaxCol; icol++)
	327	{
	328	fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
	329	fCellCount[idet][ism][irow][icol],
	330	fCellIso[idet][ism][irow][icol]);
	331	}
	332	}
	333	}
	334	}
7c11ad1f	335
1da50520	336	fclose(fpw);
7c11ad1f	337
1da50520	338	}
	339
	340	// ------------------------------------------------------------------------ //
	341
338c285a	342	Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader rawReader, TObjArray pmdddlcont)
4fc06b3a	343	{
4fc06b3a	344	// Calculates the ADC of isolated cell
7c11ad1f	345
4fc06b3a	346	const Int_t kDDL = AliDAQ::NumberOfDdls("PMD");
4fc06b3a	347	const Int_t kCellNeighbour = 6;
4fc06b3a	348	Int_t neibx[6] = {1,0,-1,-1,0,1};
	349	Int_t neiby[6] = {0,1,1,0,-1,-1};
	350
be8b7039	351	Int_t id1 = 0,jd1 = 0; //neighbour row/col
be8b7039	352	Int_t isocount = 0; //number of neighbours with 0 signal
4fc06b3a	353
0ab3a530	354	Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
7c11ad1f	355
4fc06b3a	356	for(Int_t idet = 0; idet < kDet; idet++)
	357	{
	358	for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	359	{
7c11ad1f	360	for(Int_t irow = 0; irow < kMaxRow; irow++)
4fc06b3a	361	{
	362	for(Int_t icol = 0; icol < kMaxCol; icol++)
	363	{
7c11ad1f	364	d1[idet][ismn][irow][icol] = 0.;
4fc06b3a	365	}
	366	}
	367	}
	368	}
7c11ad1f	369
4fc06b3a	370	AliPMDRawStream rawStream(rawReader);
4fc06b3a	371
1da50520	372	Int_t iddl = -1;
7c11ad1f	373
722ccc67	374	Int_t numberofDDLs = 0;
722ccc67	375
7c11ad1f	376	while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
	377	numberofDDLs++;
	378
	379	Int_t ientries = pmdddlcont->GetEntries();
	380
	381	for (Int_t ient = 0; ient < ientries; ient++)
4fc06b3a	382	{
7c11ad1f	383	AliPMDddldata pmdddl = (AliPMDddldata)pmdddlcont->UncheckedAt(ient);
	384
	385	Int_t idet = pmdddl->GetDetector();
	386	Int_t ismn = pmdddl->GetSMN();
	387	Int_t mcm = pmdddl->GetMCM();
	388	//Int_t ichno = pmdddl->GetChannel();
	389	Int_t irow = pmdddl->GetRow();
	390	Int_t icol = pmdddl->GetColumn();
	391	Int_t isig = pmdddl->GetSignal();
	392
	393	// This is the protection not to crash the code
	394
	395	if(mcm == 0) continue;
	396	if (irow < 0 \|\| icol < 0 \|\| irow > 47 \|\| icol > 95) continue;
	397
	398	// Pedestal subtraction
	399	if(fHotFlag[idet][ismn][irow][icol] == 1.0) isig = 0;
	400
	401	if (isig>0)
	402	{
	403	d1[idet][ismn][irow][icol] =
	404	(Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
	405
	406	//printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
	407
	408	fNhitCell[idet][ismn][irow][icol]++; //cell hit frequency
	409	}
	410
1da50520	411
4fc06b3a	412	}
7c11ad1f	413	pmdddlcont->Delete();
7c11ad1f	414	}
4fc06b3a	415
4fc06b3a	416	for(Int_t idet=0; idet < kDet; idet++)
7c11ad1f	417	{
4fc06b3a	418	for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
7c11ad1f	419	{
4fc06b3a	420	for(Int_t irow = 0; irow < kMaxRow; irow++)
7c11ad1f	421	{
4fc06b3a	422	for(Int_t icol = 0; icol < kMaxCol; icol++)
7c11ad1f	423	{
4fc06b3a	424	if(d1[idet][ismn][irow][icol] > 0)
7c11ad1f	425	{
4fc06b3a	426	isocount = 0;
4fc06b3a	427	for(Int_t ii = 0; ii < kCellNeighbour; ii++)
7c11ad1f	428	{
4fc06b3a	429	id1 = irow + neibx[ii];
4fc06b3a	430	jd1 = icol + neiby[ii];
d55d21ef	431	if (id1 < 0) id1 = 0;
	432	if (id1 > kMaxRow-1) id1 = kMaxRow - 1;
	433	if (jd1 < 0) jd1 = 0;
	434	if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1;
4fc06b3a	435	if(d1[idet][ismn][id1][jd1] == 0)
7c11ad1f	436	{
4fc06b3a	437	isocount++;
4fc06b3a	438	if(isocount == kCellNeighbour)
7c11ad1f	439	{
7c11ad1f	440	//fDetIso[idet] += d1[idet][ismn][irow][icol];
0ab3a530	441	fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
0ab3a530	442	fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
7c11ad1f	443	//fDetCount[idet]++;
0ab3a530	444	fSMCount[idet][ismn]++;
	445	fCellCount[idet][ismn][irow][icol]++;
	446
7c11ad1f	447	}
	448	}
	449	} // neigh cell cond.
	450	} // d>0 cond.
	451	}
	452	}
	453	}
	454	}
	455
	456	for(Int_t idet=0; idet < kDet; idet++)
	457	{
	458	for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
	459	{
	460	for(Int_t irow = 0; irow < kMaxRow; irow++)
	461	{
	462	for(Int_t icol = 0; icol < kMaxCol; icol++)
	463	{
	464	if(fNhitCell[idet][ismn][irow][icol]>0)
	465	{
	466	fCountSm[idet][ismn] += 1;
	467	fTempnhit[idet][ismn] += fNhitCell[idet][ismn][irow][icol];
	468	fTempnhitSq[idet][ismn] += fNhitCell[idet][ismn][irow][icol]
	469	*fNhitCell[idet][ismn][irow][icol];
	470	}
	471	}
	472	}
	473	}
	474	}
	475
	476	if (numberofDDLs < kDDL) return kFALSE;
722ccc67	477	return kTRUE;
7c11ad1f	478
4fc06b3a	479	}
4fc06b3a	480	// ------------------------------------------------------------------------ //
9d1486e3	481	void AliPMDCalibGain::Analyse(TTree gaintree, TTree meantree)
4fc06b3a	482	{
7c11ad1f	483	// Calculates the mean
be8b7039	484	Int_t det = 0, sm = 0, row = 0, col = 0;
be8b7039	485	Float_t gain = 0.;
7c11ad1f	486	Float_t cellmean = 0.;
7c11ad1f	487
9d1486e3	488	Float_t modmean[2][24];
	489
	490	for (Int_t idet=0; idet < 2; idet++)
	491	{
	492	for (Int_t ism = 0; ism < 24; ism++)
	493	{
	494	modmean[idet][ism] = 0.;
	495	}
	496	}
	497
7c11ad1f	498	gaintree->Branch("det",&det,"det/I");
	499	gaintree->Branch("sm",&sm,"sm/I");
	500	gaintree->Branch("row",&row,"row/I");
	501	gaintree->Branch("col",&col,"col/I");
	502	gaintree->Branch("gain",&gain,"gain/F");
	503
	504	for(Int_t idet = 0; idet < kDet; idet++)
	505	{
	506	for(Int_t ism = 0; ism < kMaxSMN; ism++)
	507	{
	508	if (fSMCount[idet][ism] > 0)
9d1486e3	509	modmean[idet][ism] = fSMIso[idet][ism]/fSMCount[idet][ism];
7c11ad1f	510	for(Int_t irow = 0; irow < kMaxRow; irow++)
	511	{
	512	for(Int_t icol = 0; icol < kMaxCol; icol++)
	513	{
	514	if (fCellCount[idet][ism][irow][icol] > 0.)
	515	{
	516	cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
	517	}
	518	det = idet;
	519	sm = ism;
	520	row = irow;
	521	col = icol;
	522	if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
	523	{
9d1486e3	524	gain = cellmean/modmean[idet][ism];
7c11ad1f	525	}
	526	else
	527	{
	528	gain = 0.;
	529	}
	530	//if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
	531	gaintree->Fill();
	532	}
	533	}
	534	}
	535	}
9d1486e3	536
	537	Float_t smmean;
	538
	539	// Writing each module mean value
	540	meantree->Branch("det",&det,"det/I");
	541	meantree->Branch("sm",&sm,"sm/I");
	542	meantree->Branch("smmean",&smmean,"row/F");
7c11ad1f	543
9d1486e3	544	for(Int_t idet = 0; idet < kDet; idet++)
	545	{
	546	for (Int_t ism = 0; ism < kMaxSMN; ism++)
	547	{
	548	det = idet;
	549	sm = ism;
	550	smmean = modmean[idet][ism];
	551	meantree->Fill();
	552	}
	553	}
	554
4fc06b3a	555	}
4fc06b3a	556	// ------------------------------------------------------------------------ //
9d1486e3	557	void AliPMDCalibGain::FindHotCell(TTree *hottree, Float_t xvar)
7c11ad1f	558	{
7c11ad1f	559	// Calculates the mean
be8b7039	560	Int_t det = 0, sm = 0, row = 0, col = 0;
	561	Float_t flag = 0.;
	562	Float_t meannhit = 0.;
	563	Float_t meanSqnhit = 0.;
	564	Float_t sigmanhit = 0.,nhitcut = 0.;
9d1486e3	565
	566	//Float_t xvar = 5.;
	567
7c11ad1f	568	hottree->Branch("det",&det,"det/I");
	569	hottree->Branch("sm",&sm,"sm/I");
	570	hottree->Branch("row",&row,"row/I");
	571	hottree->Branch("col",&col,"col/I");
	572	hottree->Branch("flag",&flag,"flag/F");
	573
	574	for(Int_t idet = 0; idet < kDet; idet++)
	575	{
	576	for(Int_t ism = 0; ism < kMaxSMN; ism++)
	577	{
	578	if (fCountSm[idet][ism]> 0)
	579	{
	580	meannhit = fTempnhit[idet][ism]/fCountSm[idet][ism];
	581	meanSqnhit = fTempnhitSq[idet][ism]/fCountSm[idet][ism];
	582	sigmanhit = sqrt(meanSqnhit-(meannhit*meannhit));
9d1486e3	583	nhitcut = meannhit + xvar*sigmanhit;
7c11ad1f	584
	585	for(Int_t irow = 0; irow < kMaxRow; irow++)
	586	{
	587	for(Int_t icol = 0; icol < kMaxCol; icol++)
	588	{
	589	det = idet;
	590	sm = ism;
	591	row = irow;
	592	col = icol;
	593
	594	if(fNhitCell[idet][ism][irow][icol] > nhitcut)
	595	{
	596	flag = 1.0;
	597	}
	598	else
	599	{
	600	flag = 0.;
	601	}
	602	hottree->Fill();
	603	}
	604
	605	}
	606	}
	607	}
	608	}
	609	}
	610