[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 "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, sm, row, col;
  Float_t mean, rms;
  
  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, sm, row, col;
  Float_t flag;

  TFile *hotmapfile = new TFile(rootFile);

  if(!hotmapfile)
    {
      printf(" NO HOTCHANNEL MAP FOUND (PMD_HOT.root) FILE IS FOUND \n");
      fHotFlag[kDet][kMaxSMN][kMaxRow][kMaxCol] = 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
  
  fpw = fopen(tempFile,"r");
  
  Float_t smcount, smiso;
  Float_t cellcount, celliso;


  for (Int_t idet = 0; idet < kDet; idet++)
    {
      for (Int_t ism = 0; ism < kMaxSMN; ism++)
	{
	  fscanf(fpw,"%d %d %f %f",&idet,&ism,&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++)
		{
		  fscanf(fpw,"%d %d %d %d %f %f",&idet,&ism,&irow,&icol,
			 &cellcount,&celliso);
		  
		  fCellCount[idet][ism][irow][icol] = cellcount;
		  fCellIso[idet][ism][irow][icol]   = celliso;
		}
	    }
	}
    }
  
  fclose(fpw);

}
// ------------------------------------------------------------------------ //
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,jd1;  //neighbour row/col
  Int_t isocount; //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, sm, row, col;
  Float_t gain;
  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, sm, row, col;
  Float_t flag;
  Float_t meannhit;
  Float_t meanSqnhit;
  Float_t sigmanhit,nhitcut;

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