[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 "AliFMDAnaParameters.h"
#include "AliFMDParameters.h"

ClassImp(AliFMDAnalysisTaskSharing)

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