[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskSharing.cxx

 
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <iostream>
#include <TMath.h>
#include "AliFMDDebug.h"
#include "AliFMDAnalysisTaskSharing.h"
#include "AliAnalysisManager.h"
#include "AliESDFMD.h"
#include "AliFMDGeometry.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliESDVertex.h"
#include "AliMultiplicity.h"
#include "AliFMDAnaParameters.h"
#include "AliFMDParameters.h"

ClassImp(AliFMDAnalysisTaskSharing)

//_____________________________________________________________________
AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
: fDebug(0),
  fESD(0x0),
  foutputESDFMD(),
  fEnergy(0),
  fNstrips(0),
  fSharedThis(kFALSE),
  fSharedPrev(kFALSE),
  fDiagList(),
  fStandalone(kTRUE),
  fEsdVertex(0),
  fStatus(kTRUE)
{
  // Default constructor
  DefineInput (0, AliESDEvent::Class());
  DefineOutput(0, AliESDFMD::Class());
  DefineOutput(1, AliESDVertex::Class());
  DefineOutput(2, AliESDEvent::Class());
  DefineOutput(3, TList::Class());
}
//_____________________________________________________________________
AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
    AliAnalysisTask(name, "AnalysisTaskFMD"),
    fDebug(0),
    fESD(0x0),
    foutputESDFMD(),
    fEnergy(0),
    fNstrips(0),
    fSharedThis(kFALSE),
    fSharedPrev(kFALSE),
    fDiagList(),
    fStandalone(kTRUE),
    fEsdVertex(0),
    fStatus(kTRUE)
{
  fStandalone = SE;
  if(fStandalone) {
    DefineInput (0, AliESDEvent::Class());
    DefineOutput(0, AliESDFMD::Class());
    DefineOutput(1, AliESDVertex::Class());
    DefineOutput(2, AliESDEvent::Class());
    DefineOutput(3, TList::Class());
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskSharing::CreateOutputObjects()
{
  if(!foutputESDFMD)
    foutputESDFMD = new AliESDFMD();
  
  if(!fEsdVertex)
    fEsdVertex    = new AliESDVertex();
  //Diagnostics
  fDiagList.SetName("Sharing diagnostics");
  for(Int_t det = 1; det<=3; det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    
    for(Int_t iring = 0;iring<nRings; iring++) {
      Char_t ringChar = (iring == 0 ? 'I' : 'O');
      TH1F* hEdist        = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
				     Form("Edist_before_sharing_FMD%d%c", det, ringChar),
				     1000,0,25);
      TH1F* hEdist_after  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
				     Form("Edist_after_sharing_FMD%d%c", det, ringChar),
				     1000,0,25);
      
      
      TH1F* hNstripsHit    = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
				     Form("N_strips_hit_FMD%d%c",det,ringChar),
				     25,0,25);
      fDiagList.Add(hEdist);
      fDiagList.Add(hEdist_after);
      fDiagList.Add(hNstripsHit);

    }
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
{
  if(fStandalone)
    fESD = (AliESDEvent*)GetInputData(0);
}
//_____________________________________________________________________
void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
{

  AliESD* old = fESD->GetAliESDOld();
  if (old) {
    fESD->CopyFromOldESD();
  }
  
  foutputESDFMD->Clear();
  
  Double_t vertex[3];
  GetVertex(vertex);
  fEsdVertex->SetXYZ(vertex);
  
  if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0) {
  
    fStatus = kFALSE;
    return;
  }
  else
    fStatus = kTRUE;
  const AliMultiplicity* testmult = fESD->GetMultiplicity();
  
  Int_t nTrackLets = testmult->GetNumberOfTracklets();
  
  if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
  else foutputESDFMD->SetUniqueID(kFALSE);
  
  AliESDFMD* fmd = fESD->GetFMDData();
  
  if (!fmd) return;
  Int_t nHits = 0;
  for(UShort_t det=1;det<=3;det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      Char_t   ring = (ir == 0 ? 'I' : 'O');
      UShort_t nsec = (ir == 0 ? 20  : 40);
      UShort_t nstr = (ir == 0 ? 512 : 256);
      
      TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
      
      for(UShort_t sec =0; sec < nsec;  sec++) {
	fSharedThis      = kFALSE;
	fSharedPrev      = kFALSE;
	fEnergy = 0;
	fNstrips = 0;
	for(UShort_t strip = 0; strip < nstr; strip++) {
	  foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
	  Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
	  
	  if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
	  
	  Double_t eta  = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip);
	  //std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<"    "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
	  
	  hEdist->Fill(mult);
	  if(fmd->IsAngleCorrected())
	    mult = mult/TMath::Cos(Eta2Theta(eta));
	  Float_t Eprev = 0;
	  Float_t Enext = 0;
	  if(strip != 0)
	    if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
	      Eprev = fmd->Multiplicity(det,ring,sec,strip-1);
	      if(fmd->IsAngleCorrected())
		Eprev = Eprev/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
	    }
	  if(strip != nstr - 1)
	    if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
	      Enext = fmd->Multiplicity(det,ring,sec,strip+1);
	      if(fmd->IsAngleCorrected())
		Enext = Enext/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
	    }
	  
	  Float_t merged_energy = GetMultiplicityOfStrip(mult,eta,Eprev,Enext,det,ring,sec,strip);

	  if(merged_energy > 0 )
	    nHits++;
	  foutputESDFMD->SetMultiplicity(det,ring,sec,strip,merged_energy);
	  foutputESDFMD->SetEta(det,ring,sec,strip,eta);
	  
	}
      }
    }
  }
  
  if(fStandalone) {
    PostData(0, foutputESDFMD); 
    PostData(1, fEsdVertex); 
    PostData(2, fESD); 
    PostData(3, &fDiagList); 
  }
}
//_____________________________________________________________________
Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
							  Float_t eta,
							  Float_t Eprev,
							  Float_t Enext,
							  UShort_t   det,
							  Char_t  ring,
							  UShort_t sec,
							  UShort_t strip) {
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 
  Float_t merged_energy = 0;
  //Float_t nParticles = 0;
  Float_t cutLow  = 0.15;
  AliFMDParameters* recopars = AliFMDParameters::Instance();
  cutLow = (4*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
  
  
  Float_t cutHigh = pars->GetMPV(det,ring,eta);// - pars->GetSigma(det,ring,eta);
  // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
  Float_t Etotal  = mult;
  
  //if(mult > 5)
  //  std::cout<<mult<<"    "<<det<<"    "<<ring<<"   "<<sec<<"    "<<strip<<std::endl;
  //Float_t slow_particle_cut = 2*pars->GetMPV(det,ring,eta);
  
  //if(recopars->IsDead(det,ring,sec,strip))
  //  std::cout<<"dead channel"<<std::endl;
  //if(foutputESDFMD->GetUniqueID() == kTRUE) {
  // Float_t mpv = pars->GetMPV(det,ring,eta);
  
  /*
  if(foutputESDFMD->GetUniqueID() == kFALSE) {
    
    if(mult > 15)
      return 0;

    if(mult > cutLow ) {
      fEnergy = fEnergy + mult;
      fNstrips++;
    }
    
    if( (Enext < cutLow && fEnergy > 0 )   || fNstrips >2 ){
      
      //if((fEnergy*TMath::Cos(Eta2Theta(eta))) > cutPart || fNstrips > 1) {
      nParticles = 1;
      merged_energy = fEnergy*TMath::Cos(Eta2Theta(eta));
      TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
      hEdist->Fill(fEnergy);
      TH1F* hNstrips = (TH1F*)fDiagList.FindObject(Form("N_strips_hit_FMD%d%c",det,ring));
      hNstrips->Fill(fNstrips);
      // std::cout<<Form("Merged signals %f %f %f into %f , %f in det %d, ring %c, sec %d, strip %d",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,det,ring,sec,strip )<<std::endl;
      
      // }
    // else
    //std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
      
      fEnergy  = 0;
      fNstrips = 0;
    
    return merged_energy;
    }
    
    
    return 0;
    
  }
  else {*/
    //std::cout<<det<<ring<<"   "<<sec<<"    "<<strip<<"   "<<cutLow<<std::endl;
  if(fSharedThis) {
    fSharedThis      = kFALSE;
    fSharedPrev      = kTRUE;
    return 0.;
  }
  
  /*  if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
    fSharedThis      = kFALSE;
    fSharedPrev      = kFALSE;
    return 0;
    }*/
  if(mult<Enext && Enext>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
    {
      fSharedThis      = kFALSE;
      fSharedPrev      = kFALSE;
      return 0;
    }
  if(mult > 15)
    {
      std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
      fSharedThis      = kFALSE;
      fSharedPrev      = kFALSE;
      return 0;
    }
  
  if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
    Etotal += Eprev;
  }
  
  if(Enext > cutLow && Enext < cutHigh ) {
    Etotal += Enext;
    fSharedThis      = kTRUE;
  }
  TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
  hEdist->Fill(Etotal);
  
  Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
  if(Etotal > 0) {
    
    merged_energy = Etotal;
    fSharedPrev      = kTRUE;
    // if(det == 1 && ring =='I')
    // std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det )<<std::endl;
  }
    else{// if(Etotal > 0) {
      //if(det == 3 && ring =='I')
      //	std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
    fSharedThis      = kFALSE;
    fSharedPrev      = kFALSE;
  }
  // merged_energy = mult;
  
  return merged_energy; 
  //}  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskSharing::GetVertex(Double_t* vertexXYZ) 
{
  const AliESDVertex* vertex = 0;
  vertex = fESD->GetPrimaryVertex();
  if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
    vertex = fESD->GetPrimaryVertexSPD();
  if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
    vertex = fESD->GetPrimaryVertexTPC();
  if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))    
    vertex = fESD->GetVertex();
  if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
    vertex->GetXYZ(vertexXYZ);
    //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
    return;
  }
  else if (fESD->GetESDTZERO()) { 
    vertexXYZ[0] = 0;
    vertexXYZ[1] = 0;
    vertexXYZ[2] = fESD->GetT0zVertex();
    
    return;
  }
  
  return;
  
}
//_____________________________________________________________________
Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {

  Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
  
  if(eta < 0)
    theta = theta-TMath::Pi();
  
  //  std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
  return theta;
  

}
//_____________________________________________________________________
Double_t AliFMDAnalysisTaskSharing::EtaFromStrip(UShort_t det, 
						Char_t ring, 
						UShort_t sector, 
						UShort_t strip, 
						Double_t zvtx)
{
   
  AliFMDGeometry* geo = AliFMDGeometry::Instance();
  
  Double_t x,y,z;
  geo->Detector2XYZ(det,ring,sector,strip,x,y,z);
  
  Double_t r = TMath::Sqrt(x*x+y*y);
  
  Double_t z_real      = z-zvtx;
  Double_t theta       = TMath::ATan2(r,z_real);
  // std::cout<<"From EtaFromStrip "<<theta<<std::endl;
  Double_t eta         =  -1*TMath::Log(TMath::Tan(0.5*theta));
  
  // std::cout<<det<<"   "<<ring<<"   "<<sector<<"   "<<strip<<"   "<<r<<"    "<<z_real<<"   "<<theta<<"    "<<eta<<std::endl;
  
  return eta;
}
//_____________________________________________________________________
//
// EOF
//
Commit	Line	Data
3bb122c7	1
	2	#include <TROOT.h>
	3	#include <TSystem.h>
	4	#include <TInterpreter.h>
	5	#include <TChain.h>
	6	#include <TFile.h>
	7	#include <TList.h>
	8	#include <iostream>
bb8a464f	9	#include <TMath.h>
bb8a464f	10	#include "AliFMDDebug.h"
3bb122c7	11	#include "AliFMDAnalysisTaskSharing.h"
	12	#include "AliAnalysisManager.h"
	13	#include "AliESDFMD.h"
bb8a464f	14	#include "AliFMDGeometry.h"
3bb122c7	15	#include "AliMCEventHandler.h"
	16	#include "AliStack.h"
	17	#include "AliESDVertex.h"
6289b3e8	18	#include "AliMultiplicity.h"
3bb122c7	19	#include "AliFMDAnaParameters.h"
bb8a464f	20	#include "AliFMDParameters.h"
3bb122c7	21
	22	ClassImp(AliFMDAnalysisTaskSharing)
	23
	24	//_____________________________________________________________________
	25	AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
	26	: fDebug(0),
	27	fESD(0x0),
8dc7c4c2	28	foutputESDFMD(),
6289b3e8	29	fEnergy(0),
6289b3e8	30	fNstrips(0),
c78bc12b	31	fSharedThis(kFALSE),
7c3e5162	32	fSharedPrev(kFALSE),
	33	fDiagList(),
	34	fStandalone(kTRUE),
bb8a464f	35	fEsdVertex(0),
bb8a464f	36	fStatus(kTRUE)
3bb122c7	37	{
	38	// Default constructor
	39	DefineInput (0, AliESDEvent::Class());
7c3e5162	40	DefineOutput(0, AliESDFMD::Class());
	41	DefineOutput(1, AliESDVertex::Class());
	42	DefineOutput(2, AliESDEvent::Class());
	43	DefineOutput(3, TList::Class());
3bb122c7	44	}
3bb122c7	45	//_____________________________________________________________________
7c3e5162	46	AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
3bb122c7	47	AliAnalysisTask(name, "AnalysisTaskFMD"),
	48	fDebug(0),
	49	fESD(0x0),
8dc7c4c2	50	foutputESDFMD(),
bb8a464f	51	fEnergy(0),
bb8a464f	52	fNstrips(0),
c78bc12b	53	fSharedThis(kFALSE),
7c3e5162	54	fSharedPrev(kFALSE),
	55	fDiagList(),
	56	fStandalone(kTRUE),
bb8a464f	57	fEsdVertex(0),
bb8a464f	58	fStatus(kTRUE)
3bb122c7	59	{
7c3e5162	60	fStandalone = SE;
	61	if(fStandalone) {
	62	DefineInput (0, AliESDEvent::Class());
	63	DefineOutput(0, AliESDFMD::Class());
	64	DefineOutput(1, AliESDVertex::Class());
	65	DefineOutput(2, AliESDEvent::Class());
	66	DefineOutput(3, TList::Class());
	67	}
3bb122c7	68	}
	69	//_____________________________________________________________________
	70	void AliFMDAnalysisTaskSharing::CreateOutputObjects()
	71	{
7c3e5162	72	if(!foutputESDFMD)
	73	foutputESDFMD = new AliESDFMD();
	74
	75	if(!fEsdVertex)
	76	fEsdVertex = new AliESDVertex();
	77	//Diagnostics
	78	fDiagList.SetName("Sharing diagnostics");
	79	for(Int_t det = 1; det<=3; det++) {
	80	Int_t nRings = (det==1 ? 1 : 2);
	81
	82	for(Int_t iring = 0;iring<nRings; iring++) {
	83	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	84	TH1F* hEdist = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
	85	Form("Edist_before_sharing_FMD%d%c", det, ringChar),
bb8a464f	86	1000,0,25);
7c3e5162	87	TH1F* hEdist_after = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
7c3e5162	88	Form("Edist_after_sharing_FMD%d%c", det, ringChar),
bb8a464f	89	1000,0,25);
	90
	91
	92	TH1F* hNstripsHit = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
	93	Form("N_strips_hit_FMD%d%c",det,ringChar),
	94	25,0,25);
7c3e5162	95	fDiagList.Add(hEdist);
7c3e5162	96	fDiagList.Add(hEdist_after);
bb8a464f	97	fDiagList.Add(hNstripsHit);
7c3e5162	98
	99	}
	100	}
3bb122c7	101	}
	102	//_____________________________________________________________________
	103	void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t /option*/)
	104	{
7c3e5162	105	if(fStandalone)
7c3e5162	106	fESD = (AliESDEvent*)GetInputData(0);
3bb122c7	107	}
	108	//_____________________________________________________________________
	109	void AliFMDAnalysisTaskSharing::Exec(Option_t /option*/)
	110	{
bb8a464f	111
3bb122c7	112	AliESD* old = fESD->GetAliESDOld();
	113	if (old) {
	114	fESD->CopyFromOldESD();
	115	}
	116
7c3e5162	117	foutputESDFMD->Clear();
3bb122c7	118
bb8a464f	119	Double_t vertex[3];
	120	GetVertex(vertex);
	121	fEsdVertex->SetXYZ(vertex);
	122
	123	if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0) {
6289b3e8	124
bb8a464f	125	fStatus = kFALSE;
	126	return;
	127	}
	128	else
	129	fStatus = kTRUE;
6289b3e8	130	const AliMultiplicity* testmult = fESD->GetMultiplicity();
	131
	132	Int_t nTrackLets = testmult->GetNumberOfTracklets();
	133
	134	if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
	135	else foutputESDFMD->SetUniqueID(kFALSE);
bb8a464f	136
3bb122c7	137	AliESDFMD* fmd = fESD->GetFMDData();
	138
	139	if (!fmd) return;
bb8a464f	140	Int_t nHits = 0;
3bb122c7	141	for(UShort_t det=1;det<=3;det++) {
	142	Int_t nRings = (det==1 ? 1 : 2);
	143	for (UShort_t ir = 0; ir < nRings; ir++) {
	144	Char_t ring = (ir == 0 ? 'I' : 'O');
	145	UShort_t nsec = (ir == 0 ? 20 : 40);
	146	UShort_t nstr = (ir == 0 ? 512 : 256);
7c3e5162	147
	148	TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
	149
8dc7c4c2	150	for(UShort_t sec =0; sec < nsec; sec++) {
	151	fSharedThis = kFALSE;
	152	fSharedPrev = kFALSE;
bb8a464f	153	fEnergy = 0;
bb8a464f	154	fNstrips = 0;
3bb122c7	155	for(UShort_t strip = 0; strip < nstr; strip++) {
7c3e5162	156	foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
3bb122c7	157	Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
bb8a464f	158
8dc7c4c2	159	if(mult == AliESDFMD::kInvalidMult \|\| mult == 0) continue;
7c3e5162	160
cc066cb9	161	Double_t eta = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip);
bb8a464f	162	//std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<" "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
6289b3e8	163
7c3e5162	164	hEdist->Fill(mult);
bb8a464f	165	if(fmd->IsAngleCorrected())
cc066cb9	166	mult = mult/TMath::Cos(Eta2Theta(eta));
8dc823cc	167	Float_t Eprev = 0;
	168	Float_t Enext = 0;
	169	if(strip != 0)
bb8a464f	170	if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
8dc823cc	171	Eprev = fmd->Multiplicity(det,ring,sec,strip-1);
bb8a464f	172	if(fmd->IsAngleCorrected())
	173	Eprev = Eprev/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
	174	}
8dc823cc	175	if(strip != nstr - 1)
bb8a464f	176	if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
	177	Enext = fmd->Multiplicity(det,ring,sec,strip+1);
	178	if(fmd->IsAngleCorrected())
	179	Enext = Enext/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
	180	}
8dc823cc	181
bb8a464f	182	Float_t merged_energy = GetMultiplicityOfStrip(mult,eta,Eprev,Enext,det,ring,sec,strip);
6289b3e8	183
bb8a464f	184	if(merged_energy > 0 )
	185	nHits++;
	186	foutputESDFMD->SetMultiplicity(det,ring,sec,strip,merged_energy);
	187	foutputESDFMD->SetEta(det,ring,sec,strip,eta);
7c3e5162	188
3bb122c7	189	}
	190	}
	191	}
	192	}
3bb122c7	193
7c3e5162	194	if(fStandalone) {
	195	PostData(0, foutputESDFMD);
	196	PostData(1, fEsdVertex);
	197	PostData(2, fESD);
	198	PostData(3, &fDiagList);
	199	}
3bb122c7	200	}
8dc823cc	201	//_____________________________________________________________________
8dc823cc	202	Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
bb8a464f	203	Float_t eta,
8dc823cc	204	Float_t Eprev,
8dc823cc	205	Float_t Enext,
bb8a464f	206	UShort_t det,
bb8a464f	207	Char_t ring,
cc066cb9	208	UShort_t sec,
cc066cb9	209	UShort_t strip) {
8dc823cc	210	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
6289b3e8	211
bb8a464f	212	Float_t merged_energy = 0;
cc066cb9	213	//Float_t nParticles = 0;
	214	Float_t cutLow = 0.15;
	215	AliFMDParameters* recopars = AliFMDParameters::Instance();
	216	cutLow = (4recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)recopars->GetDACPerMIP());
	217
	218
	219
	220	Float_t cutHigh = pars->GetMPV(det,ring,eta);// - pars->GetSigma(det,ring,eta);
6289b3e8	221	// Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
bb8a464f	222	Float_t Etotal = mult;
bb8a464f	223
6289b3e8	224	//if(mult > 5)
6289b3e8	225	// std::cout<<mult<<" "<<det<<" "<<ring<<" "<<sec<<" "<<strip<<std::endl;
cc066cb9	226	//Float_t slow_particle_cut = 2*pars->GetMPV(det,ring,eta);
6289b3e8	227
cc066cb9	228	//if(recopars->IsDead(det,ring,sec,strip))
	229	// std::cout<<"dead channel"<<std::endl;
	230	//if(foutputESDFMD->GetUniqueID() == kTRUE) {
	231	// Float_t mpv = pars->GetMPV(det,ring,eta);
	232
	233	/*
	234	if(foutputESDFMD->GetUniqueID() == kFALSE) {
6289b3e8	235
cc066cb9	236	if(mult > 15)
	237	return 0;
	238
6289b3e8	239	if(mult > cutLow ) {
	240	fEnergy = fEnergy + mult;
	241	fNstrips++;
	242	}
bb8a464f	243
cc066cb9	244	if( (Enext < cutLow && fEnergy > 0 ) \|\| fNstrips >2 ){
	245
	246	//if((fEnergy*TMath::Cos(Eta2Theta(eta))) > cutPart \|\| fNstrips > 1) {
bb8a464f	247	nParticles = 1;
	248	merged_energy = fEnergy*TMath::Cos(Eta2Theta(eta));
	249	TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
	250	hEdist->Fill(fEnergy);
	251	TH1F* hNstrips = (TH1F*)fDiagList.FindObject(Form("N_strips_hit_FMD%d%c",det,ring));
	252	hNstrips->Fill(fNstrips);
cc066cb9	253	// std::cout<<Form("Merged signals %f %f %f into %f , %f in det %d, ring %c, sec %d, strip %d",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,det,ring,sec,strip )<<std::endl;
bb8a464f	254
6289b3e8	255	// }
bb8a464f	256	// else
bb8a464f	257	//std::cout<<Form("NO HIT for %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
cc066cb9	258
	259	fEnergy = 0;
	260	fNstrips = 0;
bb8a464f	261
bb8a464f	262	return merged_energy;
cc066cb9	263	}
	264
	265
	266	return 0;
	267
bb8a464f	268	}
cc066cb9	269	else {*/
cc066cb9	270	//std::cout<<det<<ring<<" "<<sec<<" "<<strip<<" "<<cutLow<<std::endl;
c78bc12b	271	if(fSharedThis) {
	272	fSharedThis = kFALSE;
	273	fSharedPrev = kTRUE;
8dc823cc	274	return 0.;
	275	}
	276
cc066cb9	277	/* if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
6289b3e8	278	fSharedThis = kFALSE;
	279	fSharedPrev = kFALSE;
	280	return 0;
cc066cb9	281	}*/
	282	if(mult<Enext && Enext>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
	283	{
	284	fSharedThis = kFALSE;
	285	fSharedPrev = kFALSE;
	286	return 0;
	287	}
	288	if(mult > 15)
	289	{
	290	std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
	291	fSharedThis = kFALSE;
	292	fSharedPrev = kFALSE;
	293	return 0;
	294	}
bb8a464f	295
6289b3e8	296	if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
8dc823cc	297	Etotal += Eprev;
	298	}
	299
	300	if(Enext > cutLow && Enext < cutHigh ) {
	301	Etotal += Enext;
c78bc12b	302	fSharedThis = kTRUE;
8dc823cc	303	}
6289b3e8	304	TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
6289b3e8	305	hEdist->Fill(Etotal);
bb8a464f	306
6289b3e8	307	Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
6289b3e8	308	if(Etotal > 0) {
cc066cb9	309
6289b3e8	310	merged_energy = Etotal;
c78bc12b	311	fSharedPrev = kTRUE;
cc066cb9	312	// if(det == 1 && ring =='I')
cc066cb9	313	// std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det )<<std::endl;
c78bc12b	314	}
6289b3e8	315	else{// if(Etotal > 0) {
	316	//if(det == 3 && ring =='I')
	317	// std::cout<<Form("NO HIT for %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
c78bc12b	318	fSharedThis = kFALSE;
	319	fSharedPrev = kFALSE;
	320	}
6289b3e8	321	// merged_energy = mult;
8dc823cc	322
6289b3e8	323	return merged_energy;
cc066cb9	324	//}
8dc823cc	325	}
1282ce49	326	//_____________________________________________________________________
	327	void AliFMDAnalysisTaskSharing::GetVertex(Double_t* vertexXYZ)
	328	{
	329	const AliESDVertex* vertex = 0;
	330	vertex = fESD->GetPrimaryVertex();
bb8a464f	331	if(!vertex \|\| (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
	332	vertex = fESD->GetPrimaryVertexSPD();
	333	if(!vertex \|\| (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
	334	vertex = fESD->GetPrimaryVertexTPC();
	335	if(!vertex \|\| (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
	336	vertex = fESD->GetVertex();
	337	if (vertex && (vertexXYZ[0] != 0 \|\| vertexXYZ[1] != 0 \|\| vertexXYZ[2] != 0)) {
1282ce49	338	vertex->GetXYZ(vertexXYZ);
bb8a464f	339	//std::cout<<vertex->GetName()<<" "<< vertex->GetTitle() <<" "<< vertex->GetZv()<<std::endl;
1282ce49	340	return;
	341	}
	342	else if (fESD->GetESDTZERO()) {
	343	vertexXYZ[0] = 0;
	344	vertexXYZ[1] = 0;
	345	vertexXYZ[2] = fESD->GetT0zVertex();
bb8a464f	346
1282ce49	347	return;
	348	}
	349
	350	return;
bb8a464f	351
	352	}
	353	//_____________________________________________________________________
	354	Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {
	355
	356	Float_t theta = 2TMath::ATan(TMath::Exp(-1eta));
	357
	358	if(eta < 0)
	359	theta = theta-TMath::Pi();
	360
cc066cb9	361	// std::cout<<"From eta2Theta: "<<theta<<" "<<eta<<std::endl;
bb8a464f	362	return theta;
	363
	364
	365
	366	}
	367	//_____________________________________________________________________
	368	Double_t AliFMDAnalysisTaskSharing::EtaFromStrip(UShort_t det,
	369	Char_t ring,
	370	UShort_t sector,
	371	UShort_t strip,
	372	Double_t zvtx)
	373	{
	374
	375	AliFMDGeometry* geo = AliFMDGeometry::Instance();
	376
	377	Double_t x,y,z;
	378	geo->Detector2XYZ(det,ring,sector,strip,x,y,z);
	379
	380	Double_t r = TMath::Sqrt(xx+yy);
	381
	382	Double_t z_real = z-zvtx;
	383	Double_t theta = TMath::ATan2(r,z_real);
	384	// std::cout<<"From EtaFromStrip "<<theta<<std::endl;
	385	Double_t eta = -1TMath::Log(TMath::Tan(0.5theta));
	386
	387	// std::cout<<det<<" "<<ring<<" "<<sector<<" "<<strip<<" "<<r<<" "<<z_real<<" "<<theta<<" "<<eta<<std::endl;
	388
	389	return eta;
1282ce49	390	}
3bb122c7	391	//_____________________________________________________________________
	392	//
	393	// EOF
	394	//