[u/mrichter/AliRoot.git] / TPC / AliTPCclustererKr.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.                  *
 **************************************************************************/

/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */

//-----------------------------------------------------------------
//           Implementation of the TPC cluster class
//
// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
//-----------------------------------------------------------------

#include "AliTPCclustererKr.h"
#include "AliTPCclusterKr.h"
#include <vector>
#include "TObject.h"
#include "AliPadMax.h"
#include "AliSimDigits.h"
#include "AliTPCv4.h"
#include "AliTPCParam.h"
#include "AliTPCDigitsArray.h"
#include "AliTPCvtpr.h"
#include "AliTPCClustersRow.h"
#include "TTree.h"

//used in raw data finder
#include "AliTPCROC.h"
#include "AliTPCCalPad.h"
#include "AliTPCAltroMapping.h"
#include "AliTPCcalibDB.h"
#include "AliTPCRawStream.h"
#include "AliTPCRecoParam.h"
#include "AliTPCReconstructor.h"
#include "AliRawReader.h"
#include "AliTPCCalROC.h"

ClassImp(AliTPCclustererKr)


AliTPCclustererKr::AliTPCclustererKr()
  :TObject(),
  fRawData(kFALSE),
  fRowCl(0),
  fInput(0),
  fOutput(0),
  fParam(0),
  fDigarr(0),
  fRecoParam(0),
  fIsOldRCUFormat(kFALSE)
{
//
// default constructor
//
}

AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
  :TObject(),
  fRawData(kFALSE),
  fRowCl(0),
  fInput(0),
  fOutput(0),
  fParam(0),
  fDigarr(0),
  fRecoParam(0),
  fIsOldRCUFormat(kFALSE)
{
//
// copy constructor
//
  fParam = param.fParam;
  fRecoParam=param.fRecoParam;
  fIsOldRCUFormat=param.fIsOldRCUFormat;
  fRawData=param.fRawData;
  fRowCl =param.fRowCl ;
  fInput =param.fInput ;
  fOutput=param.fOutput;
  fDigarr=param.fDigarr;
} 

AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
{
  fParam = param.fParam;
  fRecoParam=param.fRecoParam;
  fIsOldRCUFormat=param.fIsOldRCUFormat;
  fRawData=param.fRawData;
  fRowCl =param.fRowCl ;
  fInput =param.fInput ;
  fOutput=param.fOutput;
  fDigarr=param.fDigarr;
  return (*this);
}

AliTPCclustererKr::~AliTPCclustererKr()
{
  //
  // destructor
  //
}

void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
{
  if (recoParam) {
    fRecoParam = recoParam;
  }else{
    //set default parameters if not specified
    fRecoParam = AliTPCReconstructor::GetRecoParam();
    if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
  }
  return;
}


////____________________________________________________________________________
////       I/O
void AliTPCclustererKr::SetInput(TTree * tree)
{
  //
  // set input tree with digits
  //
  fInput = tree;  
  if  (!fInput->GetBranch("Segment")){
    cerr<<"AliTPCclusterKr::findClusterKr(): no proper input tree !\n";
    fInput=0;
    return;
  }
}

void AliTPCclustererKr::SetOutput(TTree * tree) 
{
  //
  // set output
  //
  fOutput= tree;  
  AliTPCClustersRow clrow;
  AliTPCClustersRow *pclrow=&clrow;  
  clrow.SetClass("AliTPCclusterKr");
  clrow.SetArray(1); // to make Clones array
  fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
}

////____________________________________________________________________________
//// with new I/O
Int_t AliTPCclustererKr::finderIO()
{
  // Krypton cluster finder for simulated events from MC

  if (!fInput) { 
    Error("Digits2Clusters", "input tree not initialised");
    return 10;
  }
  
  if (!fOutput) {
    Error("Digits2Clusters", "output tree not initialised");
    return 11;
  }

  findClusterKrIO();
  return 0;
}

Int_t AliTPCclustererKr::finderIO(AliRawReader* rawReader)
{
  // Krypton cluster finder for the TPC raw data
  //
  // fParam must be defined before

  // consider noiceROC or not

  if(rawReader)fRawData=kTRUE; //set flag to data

  if (!fOutput) {
    Error("Digits2Clusters", "output tree not initialised");
    return 11;
  }

  fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
  //   used later for memory allocation

  Bool_t isAltro=kFALSE;

  AliTPCROC * roc = AliTPCROC::Instance();
  AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
  AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
  //
  AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);

  const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
  const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
  const Int_t kNS = kNIS + kNOS;//all sectors
  Int_t zeroSup = fParam->GetZeroSup();//zero suppression parameter

  //crate TPC view
  AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
  digarr->Setup(fParam);//as usually parameters

  //
  // Loop over sectors
  //
  for(Int_t sec = 0; sec < kNS; sec++) {
    AliTPCCalROC * noiseROC;
    if(noiseTPC==0x0){
      noiseROC =new AliTPCCalROC(sec);//noise=0
    }
    else{
      noiseROC = noiseTPC->GetCalROC(sec);  // noise per given sector
    }
    Int_t nRows = 0; //number of rows in sector
    Int_t nDDLs = 0; //number of DDLs
    Int_t indexDDL = 0; //DDL index
    if (sec < kNIS) {
      nRows = fParam->GetNRowLow();
      nDDLs = 2;
      indexDDL = sec * 2;
    }else {
      nRows = fParam->GetNRowUp();
      nDDLs = 4;
      indexDDL = (sec-kNIS) * 4 + kNIS * 2;
    }

    //
    // Load the raw data for corresponding DDLs
    //
    rawReader->Reset();
    input.SetOldRCUFormat(fIsOldRCUFormat);
    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);

    if(input.Next()) {
      isAltro=kTRUE;
      // Allocate memory for rows in sector (pads(depends on row) x timebins)
      for(Int_t row = 0; row < nRows; row++) {
	digarr->CreateRow(sec,row);
      }//end loop over rows
    }
    rawReader->Reset();
    input.SetOldRCUFormat(fIsOldRCUFormat);
    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);

    //
    // Begin loop over altro data
    //
    while (input.Next()) {

      //check sector consistency
      if (input.GetSector() != sec)
	AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",sec,input.GetSector()));
      
      //check row consistency
      Short_t iRow = input.GetRow();
      if (iRow < 0 || iRow >= nRows){
	AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
		      iRow, 0, nRows -1));
	continue;
      }

      //check pad consistency
      Short_t iPad = input.GetPad();
      if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(sec,iRow))) {
	AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
		      iPad, 0, roc->GetNPads(sec,iRow) ));
	continue;
      }

      //check time consistency
      Short_t iTimeBin = input.GetTime();
      if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
	//cout<<iTimeBin<<endl;
	continue;
	AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
		      iTimeBin, 0, fRecoParam->GetLastBin() -1));
      }

      //signal
      Int_t signal = input.GetSignal();
      Short_t timeWindowLow=60, timeWindowHigh=950;
      if (signal <= zeroSup || 
	  iTimeBin < timeWindowLow ||
	  iTimeBin > timeWindowHigh 
	  ) {
	digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
	continue;
      }
      if (!noiseROC) continue;
      Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
      if (noiseOnPad>2) continue; // consider noisy pad as dead
      
      if(signal<4*noiseOnPad){
	digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
	continue;
      }

//      if(signal>12){
//	cout<<" s "<<sec<<" r "<<iRow<<" p "<<iPad<<" t "<<iTimeBin<<" dig "<<signal<<endl;      
      digarr->GetRow(sec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
//      }else digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
      

    }//end of loop over altro data
  }//end of loop over sectors
  
  SetDigArr(digarr);
  if(isAltro) findClusterKrIO();
  delete digarr;

  return 0;
}

////____________________________________________________________________________
Int_t AliTPCclustererKr::findClusterKrIO()
{
  //fParam and  fDigarr must be set to run this method

  Int_t cluster_counter=0;
  const Short_t NTotalSector=fParam->GetNSector();//number of sectors
  for(Short_t sec=0; sec<NTotalSector; sec++){
    
    //vector of maxima for each sector
    std::vector<AliPadMax*> maxima_in_sector;
    
    const Short_t Nrows=fParam->GetNRow(sec);//number of rows in sector
    for(Short_t row=0; row<Nrows; row++){
      AliSimDigits *digrow;
      if(fRawData){
	digrow = (AliSimDigits*)fDigarr->GetRow(sec,row);//real data
      }else{
	digrow = (AliSimDigits*)fDigarr->LoadRow(sec,row);//MC
      }
      if(digrow){//if pointer exist
	digrow->ExpandBuffer(); //decrunch
	const Short_t npads = digrow->GetNCols();  // number of pads
	const Short_t ntime = digrow->GetNRows(); // number of timebins
	for(Short_t np=0;np<npads;np++){
	  
	  Short_t tb_max=-1;//timebin of maximum 
	  Short_t value_maximum=-1;//value of maximum in adc
	  Short_t increase_begin=-1;//timebin when increase starts
	  Short_t sum_adc=0;//sum of adc on the pad in maximum surrounding
	  bool if_increase_begin=true;//flag - check if increasing started
	  bool if_maximum=false;//flag - check if it could be maximum

	  for(Short_t nt=0;nt<ntime;nt++){
	    Short_t adc = digrow->GetDigitFast(nt,np);
	    if(adc<6){//standard was 3
	      if(if_maximum){
		if(nt-1-increase_begin<1){//at least 2 time bins, was = 1
		  tb_max=-1;
		  value_maximum=-1;
		  increase_begin=-1;
		  sum_adc=0;
		  if_increase_begin=true;
		  if_maximum=false;
		  continue;
		}
		//insert maximum, default values and set flags
		AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
								  tb_max,
								  np,
								  row),
						       increase_begin,
						       nt-1,
						       sum_adc);
		maxima_in_sector.push_back(one_maximum);
		
		tb_max=-1;
		value_maximum=-1;
		increase_begin=-1;
		sum_adc=0;
		if_increase_begin=true;
		if_maximum=false;
	      }
	      continue;
	    }
	    
	    if(if_increase_begin){
	      if_increase_begin=false;
	      increase_begin=nt;
	    }
	    
	    if(adc>value_maximum){
	      tb_max=nt;
	      value_maximum=adc;
	      if_maximum=true;
	    }
	    sum_adc+=adc;
	    if(nt==ntime-1 && if_maximum && ntime-1-increase_begin>2){//on the edge
	      //insert maximum, default values and set flags
	      AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
								tb_max,
								np,
								row),
						     increase_begin,
						     nt-1,
						     sum_adc);
	      maxima_in_sector.push_back(one_maximum);
	      
	      tb_max=-1;
	      value_maximum=-1;
	      increase_begin=-1;
	      sum_adc=0;
	      if_increase_begin=true;
	      if_maximum=false;
	      continue;
	    }
	    
	  }//end loop over timebins
	}//end loop over pads
//      }else{
//	cout<<"Pointer does not exist!!"<<endl;
      }//end if poiner exists
    }//end loop over rows
    

    Short_t max_dig=0;
    Short_t max_sum_adc=0;
    Short_t max_nt=0;
    Short_t max_np=0;
    Short_t max_row=0;
    
    for( std::vector<AliPadMax*>::iterator mp1  = maxima_in_sector.begin();
	 mp1 != maxima_in_sector.end(); ++mp1 ) {
      
      AliTPCclusterKr *tmp=new AliTPCclusterKr();
      
      Short_t n_used_pads=1;
      Int_t cluster_value=0;
      cluster_value+=(*mp1)->GetSum();
      
      max_dig      =(*mp1)->GetAdc() ;
      max_sum_adc  =(*mp1)->GetSum() ;
      max_nt       =(*mp1)->GetTime();
      max_np       =(*mp1)->GetPad() ;
      max_row      =(*mp1)->GetRow() ;
      
      AliSimDigits *digrow_tmp;
      if(fRawData){
	digrow_tmp = (AliSimDigits*)fDigarr->GetRow(sec,(*mp1)->GetRow());
      }else{
	digrow_tmp = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp1)->GetRow());
      }

      digrow_tmp->ExpandBuffer(); //decrunch
      for(Short_t tb=(*mp1)->GetBegin(); tb<((*mp1)->GetEnd())+1; tb++){
	Short_t adc_tmp = digrow_tmp->GetDigitFast(tb,(*mp1)->GetPad());
	AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp1)->GetPad(),(*mp1)->GetRow());
	tmp->fCluster.push_back(vtpr);
      }
      
      maxima_in_sector.erase(mp1);
      mp1--;
      
      for( std::vector<AliPadMax*>::iterator mp2  = maxima_in_sector.begin();
	   mp2 != maxima_in_sector.end(); ++mp2 ) {
	
	Short_t two_pad_max=4;
	Short_t two_row_max=3;
	
	if(abs(max_row - (*mp2)->GetRow())<two_row_max && 
	   abs(max_np - (*mp2)->GetPad())<two_pad_max  &&
	   abs(max_nt - (*mp2)->GetTime())<7){//was 7
	  
	  cluster_value+=(*mp2)->GetSum();

	  n_used_pads++;
	  
	  AliSimDigits *digrow_tmp1;
	  if(fRawData){
	    digrow_tmp1 = (AliSimDigits*)fDigarr->GetRow(sec,(*mp2)->GetRow());
	  }else{
	    digrow_tmp1 = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp2)->GetRow());
	  }
	  digrow_tmp1->ExpandBuffer(); //decrunch
	  
	  for(Short_t tb=(*mp2)->GetBegin(); tb<(*mp2)->GetEnd()+1; tb++){
	    Short_t adc_tmp = digrow_tmp1->GetDigitFast(tb,(*mp2)->GetPad());
	    AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp2)->GetPad(),(*mp2)->GetRow());
	      tmp->fCluster.push_back(vtpr);
	  }
	  
	  //which one is bigger
	  if( (*mp2)->GetAdc() > max_dig ){
	    max_dig      =(*mp2)->GetAdc() ;
	    max_sum_adc  =(*mp2)->GetSum() ;
	    max_nt       =(*mp2)->GetTime();
	    max_np       =(*mp2)->GetPad() ;
	    max_row      =(*mp2)->GetRow() ;
	  } else if ( (*mp2)->GetAdc() == max_dig ){
	    if( (*mp2)->GetSum() > max_sum_adc){
	      max_dig      =(*mp2)->GetAdc() ;
	      max_sum_adc  =(*mp2)->GetSum() ;
	      max_nt       =(*mp2)->GetTime();
	      max_np       =(*mp2)->GetPad() ;
	      max_row      =(*mp2)->GetRow() ;
	    }
	  }
	  maxima_in_sector.erase(mp2);
	  mp2--;
	}
      }//inside loop
      
      //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
      if(n_used_pads==1 && cluster_value/tmp->fCluster.size()<3.6)continue;
      if(n_used_pads==2 && cluster_value/tmp->fCluster.size()<3.1)continue;
      
      tmp->SetADCcluster(cluster_value);
      tmp->SetNpads(n_used_pads);
      tmp->SetMax(AliTPCvtpr(max_dig,max_nt,max_np,max_row));
      tmp->SetSec(sec);
      tmp->SetSize();
      
      cluster_counter++;
      
      //save each cluster into file

      AliTPCClustersRow *clrow= new AliTPCClustersRow();
      fRowCl=clrow;
      clrow->SetClass("AliTPCclusterKr");
      clrow->SetArray(1);
      fOutput->GetBranch("Segment")->SetAddress(&clrow);

      Int_t tmp_cluster=0;
      TClonesArray * arr = fRowCl->GetArray();
      AliTPCclusterKr * cl = new ((*arr)[tmp_cluster]) AliTPCclusterKr(*tmp);
      
      fOutput->Fill(); 
      delete clrow;

      //end of save each cluster into file adc.root
    }//outer loop

  }//end sector for
  cout<<"Number of clusters in event: "<<cluster_counter<<endl;

  return 0;
}

////____________________________________________________________________________
Commit	Line	Data
9a1e27fe	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
	17
	18	//-----------------------------------------------------------------
	19	// Implementation of the TPC cluster class
	20	//
	21	// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
	22	//-----------------------------------------------------------------
	23
	24	#include "AliTPCclustererKr.h"
	25	#include "AliTPCclusterKr.h"
	26	#include <vector>
	27	#include "TObject.h"
	28	#include "AliPadMax.h"
	29	#include "AliSimDigits.h"
	30	#include "AliTPCv4.h"
	31	#include "AliTPCParam.h"
	32	#include "AliTPCDigitsArray.h"
	33	#include "AliTPCvtpr.h"
	34	#include "AliTPCClustersRow.h"
	35	#include "TTree.h"
	36
	37	//used in raw data finder
	38	#include "AliTPCROC.h"
	39	#include "AliTPCCalPad.h"
	40	#include "AliTPCAltroMapping.h"
	41	#include "AliTPCcalibDB.h"
	42	#include "AliTPCRawStream.h"
	43	#include "AliTPCRecoParam.h"
	44	#include "AliTPCReconstructor.h"
	45	#include "AliRawReader.h"
	46	#include "AliTPCCalROC.h"
	47
	48	ClassImp(AliTPCclustererKr)
	49
	50
	51	AliTPCclustererKr::AliTPCclustererKr()
	52	:TObject(),
	53	fRawData(kFALSE),
	54	fRowCl(0),
	55	fInput(0),
	56	fOutput(0),
	57	fParam(0),
	58	fDigarr(0),
	59	fRecoParam(0),
	60	fIsOldRCUFormat(kFALSE)
	61	{
	62	//
	63	// default constructor
	64	//
65	}
66
67	AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
68	:TObject(),
69	fRawData(kFALSE),
70	fRowCl(0),
71	fInput(0),
72	fOutput(0),
73	fParam(0),
74	fDigarr(0),
75	fRecoParam(0),
76	fIsOldRCUFormat(kFALSE)
77	{
78	//
79	// copy constructor
80	//
81	fParam = param.fParam;
82	fRecoParam=param.fRecoParam;
83	fIsOldRCUFormat=param.fIsOldRCUFormat;
84	fRawData=param.fRawData;
85	fRowCl =param.fRowCl ;
86	fInput =param.fInput ;
87	fOutput=param.fOutput;
88	fDigarr=param.fDigarr;
89	}
90
91	AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
92	{
93	fParam = param.fParam;
94	fRecoParam=param.fRecoParam;
95	fIsOldRCUFormat=param.fIsOldRCUFormat;
96	fRawData=param.fRawData;
97	fRowCl =param.fRowCl ;
98	fInput =param.fInput ;
99	fOutput=param.fOutput;
100	fDigarr=param.fDigarr;
101	return (*this);
102	}
103
104	AliTPCclustererKr::~AliTPCclustererKr()
105	{
106	//
107	// destructor
108	//
109	}
110
111	void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
112	{
113	if (recoParam) {
114	fRecoParam = recoParam;
115	}else{
116	//set default parameters if not specified
117	fRecoParam = AliTPCReconstructor::GetRecoParam();
118	if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
119	}
120	return;
121	}
122
123
124	////____________________________________________________________________________
125	//// I/O
126	void AliTPCclustererKr::SetInput(TTree * tree)
127	{
128	//
129	// set input tree with digits
130	//
131	fInput = tree;
132	if (!fInput->GetBranch("Segment")){
133	cerr<<"AliTPCclusterKr::findClusterKr(): no proper input tree !\n";
134	fInput=0;
135	return;
136	}
137	}
138
139	void AliTPCclustererKr::SetOutput(TTree * tree)
140	{
141	//
142	// set output
143	//
144	fOutput= tree;
145	AliTPCClustersRow clrow;
146	AliTPCClustersRow *pclrow=&clrow;
147	clrow.SetClass("AliTPCclusterKr");
148	clrow.SetArray(1); // to make Clones array
149	fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
150	}
151
152	////____________________________________________________________________________
153	//// with new I/O
154	Int_t AliTPCclustererKr::finderIO()
155	{
156	// Krypton cluster finder for simulated events from MC
157
158	if (!fInput) {
159	Error("Digits2Clusters", "input tree not initialised");
160	return 10;
161	}
162
163	if (!fOutput) {
164	Error("Digits2Clusters", "output tree not initialised");
165	return 11;
166	}
167
168	findClusterKrIO();
169	return 0;
170	}
171
172	Int_t AliTPCclustererKr::finderIO(AliRawReader* rawReader)
173	{
174	// Krypton cluster finder for the TPC raw data
175	//
176	// fParam must be defined before
177
178	// consider noiceROC or not
179
180	if(rawReader)fRawData=kTRUE; //set flag to data
181
182	if (!fOutput) {
183	Error("Digits2Clusters", "output tree not initialised");
184	return 11;
185	}
186
187	fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
188	// used later for memory allocation
189
190	Bool_t isAltro=kFALSE;
191
192	AliTPCROC * roc = AliTPCROC::Instance();
88f51054	193	AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
9a1e27fe	194	AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
	195	//
	196	AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
	197
	198	const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
	199	const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
	200	const Int_t kNS = kNIS + kNOS;//all sectors
	201	Int_t zeroSup = fParam->GetZeroSup();//zero suppression parameter
	202
	203	//crate TPC view
	204	AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
	205	digarr->Setup(fParam);//as usually parameters
	206
	207	//
	208	// Loop over sectors
	209	//
	210	for(Int_t sec = 0; sec < kNS; sec++) {
88f51054	211	AliTPCCalROC * noiseROC;
	212	if(noiseTPC==0x0){
	213	noiseROC =new AliTPCCalROC(sec);//noise=0
	214	}
	215	else{
	216	noiseROC = noiseTPC->GetCalROC(sec); // noise per given sector
	217	}
9a1e27fe	218	Int_t nRows = 0; //number of rows in sector
	219	Int_t nDDLs = 0; //number of DDLs
	220	Int_t indexDDL = 0; //DDL index
	221	if (sec < kNIS) {
	222	nRows = fParam->GetNRowLow();
	223	nDDLs = 2;
	224	indexDDL = sec * 2;
	225	}else {
	226	nRows = fParam->GetNRowUp();
	227	nDDLs = 4;
	228	indexDDL = (sec-kNIS) * 4 + kNIS * 2;
	229	}
	230
	231	//
	232	// Load the raw data for corresponding DDLs
	233	//
	234	rawReader->Reset();
	235	input.SetOldRCUFormat(fIsOldRCUFormat);
	236	rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
	237
	238	if(input.Next()) {
	239	isAltro=kTRUE;
	240	// Allocate memory for rows in sector (pads(depends on row) x timebins)
	241	for(Int_t row = 0; row < nRows; row++) {
	242	digarr->CreateRow(sec,row);
	243	}//end loop over rows
	244	}
88f51054	245	rawReader->Reset();
88f51054	246	input.SetOldRCUFormat(fIsOldRCUFormat);
9a1e27fe	247	rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
	248
	249	//
	250	// Begin loop over altro data
	251	//
	252	while (input.Next()) {
	253
	254	//check sector consistency
	255	if (input.GetSector() != sec)
	256	AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",sec,input.GetSector()));
	257
	258	//check row consistency
	259	Short_t iRow = input.GetRow();
	260	if (iRow < 0 \|\| iRow >= nRows){
	261	AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
	262	iRow, 0, nRows -1));
	263	continue;
	264	}
	265
	266	//check pad consistency
	267	Short_t iPad = input.GetPad();
	268	if (iPad < 0 \|\| iPad >= (Short_t)(roc->GetNPads(sec,iRow))) {
	269	AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
	270	iPad, 0, roc->GetNPads(sec,iRow) ));
	271	continue;
	272	}
	273
	274	//check time consistency
	275	Short_t iTimeBin = input.GetTime();
	276	if ( iTimeBin < fRecoParam->GetFirstBin() \|\| iTimeBin >= fRecoParam->GetLastBin()){
	277	//cout<<iTimeBin<<endl;
	278	continue;
	279	AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
	280	iTimeBin, 0, fRecoParam->GetLastBin() -1));
	281	}
	282
	283	//signal
	284	Int_t signal = input.GetSignal();
88f51054	285	Short_t timeWindowLow=60, timeWindowHigh=950;
	286	if (signal <= zeroSup \|\|
	287	iTimeBin < timeWindowLow \|\|
	288	iTimeBin > timeWindowHigh
	289	) {
9a1e27fe	290	digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
88f51054	291	continue;
9a1e27fe	292	}
88f51054	293	if (!noiseROC) continue;
	294	Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
	295	if (noiseOnPad>2) continue; // consider noisy pad as dead
	296
	297	if(signal<4*noiseOnPad){
	298	digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
	299	continue;
	300	}
	301
	302	// if(signal>12){
	303	// cout<<" s "<<sec<<" r "<<iRow<<" p "<<iPad<<" t "<<iTimeBin<<" dig "<<signal<<endl;
	304	digarr->GetRow(sec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
	305	// }else digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
	306
	307
9a1e27fe	308	}//end of loop over altro data
	309	}//end of loop over sectors
	310
	311	SetDigArr(digarr);
	312	if(isAltro) findClusterKrIO();
	313	delete digarr;
	314
	315	return 0;
	316	}
	317
	318	////____________________________________________________________________________
	319	Int_t AliTPCclustererKr::findClusterKrIO()
	320	{
	321	//fParam and fDigarr must be set to run this method
	322
	323	Int_t cluster_counter=0;
	324	const Short_t NTotalSector=fParam->GetNSector();//number of sectors
	325	for(Short_t sec=0; sec<NTotalSector; sec++){
	326
	327	//vector of maxima for each sector
	328	std::vector<AliPadMax*> maxima_in_sector;
	329
	330	const Short_t Nrows=fParam->GetNRow(sec);//number of rows in sector
	331	for(Short_t row=0; row<Nrows; row++){
	332	AliSimDigits *digrow;
	333	if(fRawData){
	334	digrow = (AliSimDigits*)fDigarr->GetRow(sec,row);//real data
	335	}else{
	336	digrow = (AliSimDigits*)fDigarr->LoadRow(sec,row);//MC
	337	}
	338	if(digrow){//if pointer exist
	339	digrow->ExpandBuffer(); //decrunch
	340	const Short_t npads = digrow->GetNCols(); // number of pads
	341	const Short_t ntime = digrow->GetNRows(); // number of timebins
	342	for(Short_t np=0;np<npads;np++){
	343
	344	Short_t tb_max=-1;//timebin of maximum
	345	Short_t value_maximum=-1;//value of maximum in adc
	346	Short_t increase_begin=-1;//timebin when increase starts
	347	Short_t sum_adc=0;//sum of adc on the pad in maximum surrounding
88f51054	348	bool if_increase_begin=true;//flag - check if increasing started
9a1e27fe	349	bool if_maximum=false;//flag - check if it could be maximum
	350
	351	for(Short_t nt=0;nt<ntime;nt++){
	352	Short_t adc = digrow->GetDigitFast(nt,np);
88f51054	353	if(adc<6){//standard was 3
9a1e27fe	354	if(if_maximum){
88f51054	355	if(nt-1-increase_begin<1){//at least 2 time bins, was = 1
9a1e27fe	356	tb_max=-1;
	357	value_maximum=-1;
	358	increase_begin=-1;
	359	sum_adc=0;
	360	if_increase_begin=true;
	361	if_maximum=false;
	362	continue;
	363	}
	364	//insert maximum, default values and set flags
	365	AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
	366	tb_max,
	367	np,
	368	row),
	369	increase_begin,
	370	nt-1,
	371	sum_adc);
	372	maxima_in_sector.push_back(one_maximum);
	373
	374	tb_max=-1;
	375	value_maximum=-1;
	376	increase_begin=-1;
	377	sum_adc=0;
	378	if_increase_begin=true;
	379	if_maximum=false;
	380	}
	381	continue;
	382	}
	383
	384	if(if_increase_begin){
	385	if_increase_begin=false;
	386	increase_begin=nt;
	387	}
	388
	389	if(adc>value_maximum){
	390	tb_max=nt;
	391	value_maximum=adc;
	392	if_maximum=true;
	393	}
	394	sum_adc+=adc;
	395	if(nt==ntime-1 && if_maximum && ntime-1-increase_begin>2){//on the edge
	396	//insert maximum, default values and set flags
	397	AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
	398	tb_max,
	399	np,
	400	row),
	401	increase_begin,
	402	nt-1,
	403	sum_adc);
	404	maxima_in_sector.push_back(one_maximum);
	405
	406	tb_max=-1;
	407	value_maximum=-1;
	408	increase_begin=-1;
	409	sum_adc=0;
	410	if_increase_begin=true;
	411	if_maximum=false;
	412	continue;
	413	}
	414
	415	}//end loop over timebins
	416	}//end loop over pads
	417	// }else{
	418	// cout<<"Pointer does not exist!!"<<endl;
	419	}//end if poiner exists
420	}//end loop over rows
421
422
423	Short_t max_dig=0;
424	Short_t max_sum_adc=0;
425	Short_t max_nt=0;
426	Short_t max_np=0;
427	Short_t max_row=0;
428
429	for( std::vector<AliPadMax*>::iterator mp1 = maxima_in_sector.begin();
430	mp1 != maxima_in_sector.end(); ++mp1 ) {
431
432	AliTPCclusterKr *tmp=new AliTPCclusterKr();
433
434	Short_t n_used_pads=1;
88f51054	435	Int_t cluster_value=0;
9a1e27fe	436	cluster_value+=(*mp1)->GetSum();
	437
	438	max_dig =(*mp1)->GetAdc() ;
	439	max_sum_adc =(*mp1)->GetSum() ;
	440	max_nt =(*mp1)->GetTime();
	441	max_np =(*mp1)->GetPad() ;
	442	max_row =(*mp1)->GetRow() ;
	443
	444	AliSimDigits *digrow_tmp;
	445	if(fRawData){
	446	digrow_tmp = (AliSimDigits)fDigarr->GetRow(sec,(mp1)->GetRow());
	447	}else{
	448	digrow_tmp = (AliSimDigits)fDigarr->LoadRow(sec,(mp1)->GetRow());
	449	}
	450
	451	digrow_tmp->ExpandBuffer(); //decrunch
	452	for(Short_t tb=(mp1)->GetBegin(); tb<((mp1)->GetEnd())+1; tb++){
	453	Short_t adc_tmp = digrow_tmp->GetDigitFast(tb,(*mp1)->GetPad());
	454	AliTPCvtpr vtpr=new AliTPCvtpr(adc_tmp,tb,(mp1)->GetPad(),(*mp1)->GetRow());
	455	tmp->fCluster.push_back(vtpr);
	456	}
	457
	458	maxima_in_sector.erase(mp1);
	459	mp1--;
	460
	461	for( std::vector<AliPadMax*>::iterator mp2 = maxima_in_sector.begin();
	462	mp2 != maxima_in_sector.end(); ++mp2 ) {
	463
	464	Short_t two_pad_max=4;
	465	Short_t two_row_max=3;
	466
	467	if(abs(max_row - (*mp2)->GetRow())<two_row_max &&
	468	abs(max_np - (*mp2)->GetPad())<two_pad_max &&
88f51054	469	abs(max_nt - (*mp2)->GetTime())<7){//was 7
9a1e27fe	470
9a1e27fe	471	cluster_value+=(*mp2)->GetSum();
88f51054	472
9a1e27fe	473	n_used_pads++;
	474
	475	AliSimDigits *digrow_tmp1;
	476	if(fRawData){
	477	digrow_tmp1 = (AliSimDigits)fDigarr->GetRow(sec,(mp2)->GetRow());
	478	}else{
	479	digrow_tmp1 = (AliSimDigits)fDigarr->LoadRow(sec,(mp2)->GetRow());
	480	}
	481	digrow_tmp1->ExpandBuffer(); //decrunch
	482
	483	for(Short_t tb=(mp2)->GetBegin(); tb<(mp2)->GetEnd()+1; tb++){
	484	Short_t adc_tmp = digrow_tmp1->GetDigitFast(tb,(*mp2)->GetPad());
	485	AliTPCvtpr vtpr=new AliTPCvtpr(adc_tmp,tb,(mp2)->GetPad(),(*mp2)->GetRow());
	486	tmp->fCluster.push_back(vtpr);
	487	}
	488
	489	//which one is bigger
	490	if( (*mp2)->GetAdc() > max_dig ){
	491	max_dig =(*mp2)->GetAdc() ;
	492	max_sum_adc =(*mp2)->GetSum() ;
	493	max_nt =(*mp2)->GetTime();
	494	max_np =(*mp2)->GetPad() ;
	495	max_row =(*mp2)->GetRow() ;
	496	} else if ( (*mp2)->GetAdc() == max_dig ){
	497	if( (*mp2)->GetSum() > max_sum_adc){
	498	max_dig =(*mp2)->GetAdc() ;
	499	max_sum_adc =(*mp2)->GetSum() ;
	500	max_nt =(*mp2)->GetTime();
	501	max_np =(*mp2)->GetPad() ;
	502	max_row =(*mp2)->GetRow() ;
	503	}
	504	}
	505	maxima_in_sector.erase(mp2);
	506	mp2--;
	507	}
	508	}//inside loop
	509
	510	//through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
	511	if(n_used_pads==1 && cluster_value/tmp->fCluster.size()<3.6)continue;
	512	if(n_used_pads==2 && cluster_value/tmp->fCluster.size()<3.1)continue;
	513
	514	tmp->SetADCcluster(cluster_value);
	515	tmp->SetNpads(n_used_pads);
	516	tmp->SetMax(AliTPCvtpr(max_dig,max_nt,max_np,max_row));
	517	tmp->SetSec(sec);
	518	tmp->SetSize();
	519
	520	cluster_counter++;
	521
	522	//save each cluster into file
	523
	524	AliTPCClustersRow *clrow= new AliTPCClustersRow();
	525	fRowCl=clrow;
	526	clrow->SetClass("AliTPCclusterKr");
	527	clrow->SetArray(1);
	528	fOutput->GetBranch("Segment")->SetAddress(&clrow);
	529
	530	Int_t tmp_cluster=0;
	531	TClonesArray * arr = fRowCl->GetArray();
	532	AliTPCclusterKr * cl = new ((arr)[tmp_cluster]) AliTPCclusterKr(tmp);
	533
	534	fOutput->Fill();
	535	delete clrow;
	536
537	//end of save each cluster into file adc.root
538	}//outer loop
539
540	}//end sector for
541	cout<<"Number of clusters in event: "<<cluster_counter<<endl;
542
543	return 0;
544	}
545
546	////____________________________________________________________________________