[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDMCTrackDensity.cxx

#include "AliFMDMCTrackDensity.h"
#include <AliESDFMD.h>
#include <AliMCEvent.h>
#include <AliTrackReference.h>
#include <AliStack.h>
#include <TMath.h>
#include "AliFMDStripIndex.h"
#include "AliFMDFloatMap.h"
#include <AliLog.h>
#include <TH2D.h>
#include <TH1D.h>
#include <TList.h>
#include <TROOT.h>
#include <iostream>
#include "AliGenHijingEventHeader.h"
#include <TF1.h>
#include <TGraph.h>

//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity()
  : TNamed(), 
    fUseOnlyPrimary(false), 
    fMaxConsequtiveStrips(3), 
    fNr(0), 
    fNt(0),
    fBinFlow(0), 
    fEtaBinFlow(0),
    fPhiBinFlow(0),
    fDebug(false)
{
  // Default constructor 
}

//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
  : TNamed("fmdMCTrackDensity","fmdMCTrackDensity"), 
    fUseOnlyPrimary(false), 
    fMaxConsequtiveStrips(3), 
    fNr(0), 
    fNt(0),
    fBinFlow(0), 
    fEtaBinFlow(0),
    fPhiBinFlow(0),
    fDebug(false)
{
  // Normal constructor constructor 
}

//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
  : TNamed(o),
    fUseOnlyPrimary(o.fUseOnlyPrimary), 
    fMaxConsequtiveStrips(o.fMaxConsequtiveStrips), 
    fNr(o.fNr), 
    fNt(o.fNt),
    fBinFlow(o.fBinFlow), 
    fEtaBinFlow(o.fEtaBinFlow),
    fPhiBinFlow(o.fPhiBinFlow),
    fDebug(o.fDebug)
{
  // Normal constructor constructor 
}

//____________________________________________________________________
AliFMDMCTrackDensity&
AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
{
  // Assignment operator 
  if (&o == this) return *this; 
  TNamed::operator=(o);
  fUseOnlyPrimary       = o.fUseOnlyPrimary;
  fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
  fNr                   = o.fNr;
  fNt                   = o.fNt;
  fBinFlow              = o.fBinFlow;
  fEtaBinFlow           = o.fEtaBinFlow;
  fPhiBinFlow           = o.fPhiBinFlow;
  fDebug                = o.fDebug;
  return *this;
}

//____________________________________________________________________
void
AliFMDMCTrackDensity::DefineOutput(TList* l)
{
  fNr = new TH1D("clusterRefs", "# track references per cluster",
		 21, -.5, 20.5);
  fNr->SetXTitle("# of track references/cluster");
  fNr->SetDirectory(0);
  l->Add(fNr);

  fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
  fNt->SetXTitle("Cluster size [strips]");
  fNt->SetDirectory(0);
  l->Add(fNt);

  fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
		      200, -5, 5, 40, -180, 180);
  fBinFlow->SetXTitle("#Delta#eta");
  fBinFlow->SetYTitle("#Delta#varphi");
  fBinFlow->SetDirectory(0);
  l->Add(fBinFlow);

  fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta", 
			 200, -4, 6, 200, -5, 5);
  fEtaBinFlow->SetXTitle("#eta of strip");
  fEtaBinFlow->SetYTitle("#Delta#eta");
  fEtaBinFlow->SetDirectory(0);
  l->Add(fEtaBinFlow);

  fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi", 
			 40, 0, 360, 40, -180, 180);
  fPhiBinFlow->SetXTitle("#varphi of strip");
  fPhiBinFlow->SetYTitle("#Delta#varphi");
  fPhiBinFlow->SetDirectory(0);
  l->Add(fPhiBinFlow);
}

//____________________________________________________________________
void
AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
				    const AliMCParticle* mother, 
				    Int_t                longest,
				    Double_t             vz,
				    UShort_t             nC, 
				    UShort_t             nT,
				    AliESDFMD&           output,
                                    Double_t             w) const
{
  // Store a particle. 
  if (longest < 0) return;

  AliTrackReference* ref = particle->GetTrackReference(longest);
  if (!ref) return;
    
  // Get the detector coordinates 
  UShort_t d, s, t;
  Char_t r;
  AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);

  // Check if we have value already 
  Double_t old = output.Multiplicity(d,r,s,t);

  // If invalid, force it valid 
  if (old == AliESDFMD::kInvalidMult) old = 0;

  // Increment count 
  output.SetMultiplicity(d,r,s,t,old+w);

  // Get track-reference stuff 
  Double_t x      = ref->X();
  Double_t y      = ref->Y();
  Double_t z      = ref->Z()-vz;
  Double_t rr     = TMath::Sqrt(x*x+y*y);
  Double_t thetaR = TMath::ATan2(rr,z);
  Double_t phiR   = TMath::ATan2(y,x);
  Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
  
  // Correct angle and convert to degrees 
  if (thetaR < 0) thetaR += 2*TMath::Pi();
  thetaR *= 180. / TMath::Pi();
  if (phiR < 0) phiR += 2*TMath::Pi();
  phiR *= 180. / TMath::Pi();

  // Fill histograms 
  fNr->Fill(nC);
  fNt->Fill(nT);

  const AliMCParticle* mp = (mother ? mother : particle);
  Double_t dEta = mp->Eta() - etaR;
  Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
  if (dPhi >  180) dPhi -= 360;
  if (dPhi < -180) dPhi += 360;
  fBinFlow->Fill(dEta, dPhi);
  fEtaBinFlow->Fill(etaR, dEta);
  fPhiBinFlow->Fill(phiR, dPhi);

  if (fDebug) 
    Info("StoreParticle", "Store @ FMD%d%c[%2d,%3d] nStrips=%3d, "
	 "dEta=%7.4f, dPhi=%d",
	 d, r, s, t, nT, dEta, int(dPhi+.5));
}

//____________________________________________________________________
Double_t
AliFMDMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref,
				       Double_t vz) const
{
  // Get the incidient angle of the track reference. 
  Double_t x    = ref->X();
  Double_t y    = ref->Y();
  Double_t z    = ref->Z()-vz;
  Double_t rr   = TMath::Sqrt(x*x+y*y);
  Double_t theta= TMath::ATan2(rr,z);
  Double_t ang  = TMath::Abs(TMath::Pi()-theta);
  return ang;
}
				    
//____________________________________________________________________
const AliMCParticle*
AliFMDMCTrackDensity::GetMother(Int_t     iTr,
				const AliMCEvent& event) const
{
  // 
  // Track down primary mother 
  // 
  Int_t i  = iTr;
  do { 
    const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
    if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
    
    i = p->GetMother();
  } while (i > 0);

  return 0;
}  

// #define USE_FLOW_WEIGHTS 1
//____________________________________________________________________
Bool_t
AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input, 
				const AliMCEvent& event,
				Double_t          vz,
				AliESDFMD&        output, 
				TH2D*             primary)
{
  // 
  // Filter the input kinematics and track references, using 
  // some of the ESD information
  // 
  // Parameters:
  //    input   Input ESD event
  //    event   Input MC event
  //    vz      Vertex position 
  //    output  Output ESD-like object
  //    primary Per-event histogram of primaries 
  //
  // Return:
  //    True on succes, false otherwise 
  //
  output.Clear();

  // Copy eta values to output 
  for (UShort_t ed = 1; ed <= 3; ed++) { 
    UShort_t nq = (ed == 1 ? 1 : 2);
    for (UShort_t eq = 0; eq < nq; eq++) {
      Char_t   er = (eq == 0 ? 'I' : 'O');
      UShort_t ns = (eq == 0 ?  20 :  40);
      UShort_t nt = (eq == 0 ? 512 : 256);
      for (UShort_t es = 0; es < ns; es++) 
	for (UShort_t et = 0; et < nt; et++) 
	  output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
    }
  }

#ifdef USE_FLOW_WEIGHTS
  AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
                (event.GenEventHeader());
  Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
  Double_t b = (hd ? hd->ImpactParameter() : -1 );
#endif

  AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
  Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
  Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
  for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
    AliMCParticle* particle = 
      static_cast<AliMCParticle*>(event.GetTrack(iTr));
    
    // Check the returned particle 
    if (!particle) continue;
    
    // Check if this charged and a primary 
    Bool_t isCharged = particle->Charge() != 0;
    if (!isCharged) continue;
    Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);

    // Pseudo rapidity and azimuthal angle 
    Double_t eta = particle->Eta();
    Double_t phi = particle->Phi();
    
    // Fill 'dn/deta' histogram 
    if (isPrimary && iTr < nPrim) {
      if (primary) primary->Fill(eta, phi);
    }

    // Bail out if we're only processing primaries - perhaps we should
    // track back to the original primary?
    if (fUseOnlyPrimary && !isPrimary) continue;

    Int_t    nTrRef   = particle->GetNumberOfTrackReferences();
    Int_t    longest  = -1;
    Double_t angle    = 0;
    UShort_t oD       = 0, oS = 1024, oT = 1024, ooT=1024;
    Char_t   oR       = '\0';
    UShort_t nC       = 0;
    UShort_t nT       = 0;
    // Double_t oTheta= 0;
    for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
      AliTrackReference* ref = particle->GetTrackReference(iTrRef);
      
      // Check existence 
      if (!ref) continue;

      // Check that we hit an FMD element 
      if (ref->DetectorId() != AliTrackReference::kFMD) 
	continue;

      // Get the detector coordinates 
      UShort_t d, s, t;
      Char_t r;
      AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);

      // Calculate length of last and second to last strips. 

      // IF base of cluster not set, set it here. 
      if (ooT == 1024) ooT = t;

      // Calculate distance of previous reference to base of cluster 
      nT = TMath::Abs(t - ooT) + 1;

      // Count number of track refs in this sector 
      nC++;

      Bool_t used = false;
      // If this is a new detector/ring, then reset the other one 
      // Check if we have a valid old detectorm ring, and sector 
      if (oD >  0 && oR != '\0' && oS != 1024) {
	// New detector, new ring, or new sector 
	if (d != oD   || r != oR   || s != oS) {
	  if (fDebug) Info("Process", "New because new sector");
	  used = true;
	}
      }
      if (nT > fMaxConsequtiveStrips) {
	if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)", 
			 nT, t, oT, ooT);
	used = true;
      }
      if (used) {
	if (fDebug) 
	  Info("Process", "I=%3d L=%3d D=%d (was %d), R=%c (was %c), "
	       "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
	       iTr, longest, d, oD, r, oR, s, oS, t, oT, nT, 
	       fMaxConsequtiveStrips);
	Int_t nnT   = TMath::Abs(oT - ooT) + 1;
	const AliMCParticle* mother = GetMother(iTr, event);
	
#ifdef USE_FLOW_WEIGHTS
        phi = (mother ? mother->Phi() : particle->Phi());
        eta = (mother ? mother->Eta() : particle->Eta());
        Double_t pt = (mother ? mother->Pt() : particle->Pt());
        Int_t   id = (mother ? mother->PdgCode() : 2212);
	Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
#else
        Double_t weight = 1; // // 
#endif
        StoreParticle(particle, mother, longest, vz, nC, nnT, output, weight);
	longest = -1;
	angle   = 0;
	nC  = 1;    // Reset track-ref counter - we have this->1
	nT  = 0;    // Reset cluster size 
	oD  = 0;    // Reset detector
	oR  = '\0'; // Reset ring 
	oS  = 1024; // Reset sector 
	oT  = 1024; // Reset old strip
	ooT = t;    // Reset base 
      }
      else if (fDebug) {
	if (ooT == t) 
	  Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]", 
	       d, r, s, t);
	else 
	  Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
	       "length=%3d (now in %3d, previous %3d)", 
	       d, r, s, ooT, nT, t, oT);
      }
	
      oD = d;
      oR = r;
      oS = s;
      oT = t;
      nT = TMath::Abs(t-ooT)+1;

      // The longest passage is determined through the angle 
      Double_t ang  = GetTrackRefTheta(ref, vz);
      if (ang > angle) {
	longest = iTrRef;
	angle   = ang;
      }
      // oTheta = ang;
    } // Loop over track references
    if (longest < 0) continue; // Nothing found

    // Get the reference corresponding to the longest path through the detector
    if (fDebug) 
      Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)", 
	   iTr, longest, nT, fMaxConsequtiveStrips);
    const AliMCParticle* mother = GetMother(iTr, event);

#ifdef USE_FLOW_WEIGHTS
    phi = (mother ? mother->Phi() : particle->Phi());
    eta = (mother ? mother->Eta() : particle->Eta());
    Double_t pt = (mother ? mother->Pt() : particle->Pt());
    Int_t    id = (mother ? mother->PdgCode() : 2212);
	Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
#else
        Double_t weight = 1; // // 
#endif
    StoreParticle(particle, mother, longest, vz, nC, nT, output, weight);
  } // Loop over tracks
  return kTRUE;
}

//____________________________________________________________________
Double_t
AliFMDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b, 
				      Double_t phi, Double_t rp, Int_t id) const
{
  static TF1 gaus = TF1("gaus", "gaus", -6, 6);
  gaus.SetParameters(0.1, 0., 9);
  //  gaus.SetParameters(0.1, 0., 3);
  //  gaus.SetParameters(0.1, 0., 15);
  
  const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350,
					    0.45, 0.55, 0.650, 
					    0.75, 0.85, 0.950,
					    1.10, 1.30, 1.500,
					    1.70, 1.90, 2.250,
					    2.75, 3.25, 3.750,
					    4.50};
  const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400,
					    0.059911,0.073516,
					    0.089756,0.105486,
					    0.117391,0.128199,
					    0.138013,0.158271,
					    0.177726,0.196383,
					    0.208277,0.216648,
					    0.242954,0.249961,
					    0.240131,0.269006,
					    0.207796};
  const Int_t nPointsCumulant2nd4050ALICE = 
    sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
  static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
#if 0
  const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35,
					    0.45,0.550,0.65,
					    0.75,0.850,0.95,
					    1.10,1.300,1.50,
					    1.70,1.900,2.25,
					    2.75,3.250,3.75,
					    4.50,5.500,7.00,
					    9.000000};
  const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,
					    0.048566,0.061083,
					    0.070910,0.078831,
					    0.091396,0.102026,
					    0.109691,0.124449,
					    0.139819,0.155561,
					    0.165701,0.173678,
					    0.191149,0.202015,
					    0.204540,0.212560,
					    0.195885,0.000000,
					    0.000000,0.000000};
#endif
  const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350,
					    0.45,0.55,0.650,
					    0.75,0.85,0.950,
					    1.10,1.30,1.500,
					    1.70,1.90,2.250,
					    2.75,3.25,3.750,
					    4.50};
  const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,
					    0.049824,0.066662,
					    0.075856,0.081583,
					    0.099778,0.104674,
					    0.118545,0.131874,
					    0.152959,0.155348,
					    0.169751,0.179052,
					    0.178532,0.198851,
					    0.185737,0.239901,
					    0.186098};
  const Int_t nPointsCumulant4th4050ALICE = 
    sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
  static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE, 
			       xCumulant4th4050ALICE, 
			       yCumulant4th4050ALICE);

  const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75,
						       4.225,
						       5.965,
						       7.765,
						       9.215,
						       10.46,
						       11.565,
						       12.575};
  const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,
						       0.055818,0.073137,
						       0.083898,0.086690,
						       0.082040,0.077777};
  const Int_t nPointsCumulant4thTPCrefMultTPConlyAll = 
    sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
  TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll,
		   xCumulant4thTPCrefMultTPConlyAll,
		   yCumulant4thTPCrefMultTPConlyAll);


  Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 + 
					     alicePointsPt4.Eval(pt) * 0.5) 
		     * (aliceCent.Eval(b) / aliceCent.Eval(10.46)) 
		     * 2. * TMath::Cos(2. * (phi - rp)));
  if      (TMath::Abs(id) == 211)  weight *= 1.3; //pion flow
  else if (TMath::Abs(id) == 2212) weight *= 1.0;  //proton flow
  else                             weight *= 0.7;
  
  return weight;
}
//____________________________________________________________________
void
AliFMDMCTrackDensity::Print(Option_t* /*option*/) const 
{
  char ind[gROOT->GetDirLevel()+1];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << ClassName() << ": " << GetName() << '\n'
	    << std::boolalpha 
	    << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
	    << ind << " Max cluster size:       " << fMaxConsequtiveStrips
	    << std::noboolalpha << std::endl;
  
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
f68d9069	1	#include "AliFMDMCTrackDensity.h"
	2	#include <AliESDFMD.h>
	3	#include <AliMCEvent.h>
	4	#include <AliTrackReference.h>
	5	#include <AliStack.h>
	6	#include <TMath.h>
	7	#include "AliFMDStripIndex.h"
	8	#include "AliFMDFloatMap.h"
	9	#include <AliLog.h>
	10	#include <TH2D.h>
	11	#include <TH1D.h>
	12	#include <TList.h>
	13	#include <TROOT.h>
	14	#include <iostream>
2b556440	15	#include "AliGenHijingEventHeader.h"
	16	#include <TF1.h>
	17	#include <TGraph.h>
f68d9069	18
	19	//____________________________________________________________________
	20	AliFMDMCTrackDensity::AliFMDMCTrackDensity()
	21	: TNamed(),
	22	fUseOnlyPrimary(false),
	23	fMaxConsequtiveStrips(3),
	24	fNr(0),
	25	fNt(0),
	26	fBinFlow(0),
	27	fEtaBinFlow(0),
	28	fPhiBinFlow(0),
	29	fDebug(false)
	30	{
	31	// Default constructor
	32	}
	33
	34	//____________________________________________________________________
	35	AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
	36	: TNamed("fmdMCTrackDensity","fmdMCTrackDensity"),
	37	fUseOnlyPrimary(false),
	38	fMaxConsequtiveStrips(3),
	39	fNr(0),
	40	fNt(0),
	41	fBinFlow(0),
	42	fEtaBinFlow(0),
	43	fPhiBinFlow(0),
	44	fDebug(false)
	45	{
	46	// Normal constructor constructor
	47	}
	48
	49	//____________________________________________________________________
	50	AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
	51	: TNamed(o),
	52	fUseOnlyPrimary(o.fUseOnlyPrimary),
	53	fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
	54	fNr(o.fNr),
	55	fNt(o.fNt),
	56	fBinFlow(o.fBinFlow),
	57	fEtaBinFlow(o.fEtaBinFlow),
	58	fPhiBinFlow(o.fPhiBinFlow),
	59	fDebug(o.fDebug)
	60	{
	61	// Normal constructor constructor
	62	}
	63
	64	//____________________________________________________________________
	65	AliFMDMCTrackDensity&
	66	AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
	67	{
	68	// Assignment operator
d015ecfe	69	if (&o == this) return *this;
f68d9069	70	TNamed::operator=(o);
	71	fUseOnlyPrimary = o.fUseOnlyPrimary;
	72	fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
	73	fNr = o.fNr;
	74	fNt = o.fNt;
	75	fBinFlow = o.fBinFlow;
	76	fEtaBinFlow = o.fEtaBinFlow;
	77	fPhiBinFlow = o.fPhiBinFlow;
	78	fDebug = o.fDebug;
	79	return *this;
	80	}
	81
	82	//____________________________________________________________________
	83	void
	84	AliFMDMCTrackDensity::DefineOutput(TList* l)
	85	{
	86	fNr = new TH1D("clusterRefs", "# track references per cluster",
	87	21, -.5, 20.5);
	88	fNr->SetXTitle("# of track references/cluster");
	89	fNr->SetDirectory(0);
	90	l->Add(fNr);
	91
	92	fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
	93	fNt->SetXTitle("Cluster size [strips]");
	94	fNt->SetDirectory(0);
	95	l->Add(fNt);
	96
	97	fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
	98	200, -5, 5, 40, -180, 180);
	99	fBinFlow->SetXTitle("#Delta#eta");
	100	fBinFlow->SetYTitle("#Delta#varphi");
	101	fBinFlow->SetDirectory(0);
	102	l->Add(fBinFlow);
	103
	104	fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta",
	105	200, -4, 6, 200, -5, 5);
	106	fEtaBinFlow->SetXTitle("#eta of strip");
	107	fEtaBinFlow->SetYTitle("#Delta#eta");
	108	fEtaBinFlow->SetDirectory(0);
	109	l->Add(fEtaBinFlow);
	110
	111	fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi",
	112	40, 0, 360, 40, -180, 180);
	113	fPhiBinFlow->SetXTitle("#varphi of strip");
	114	fPhiBinFlow->SetYTitle("#Delta#varphi");
	115	fPhiBinFlow->SetDirectory(0);
	116	l->Add(fPhiBinFlow);
	117	}
	118
	119	//____________________________________________________________________
	120	void
	121	AliFMDMCTrackDensity::StoreParticle(AliMCParticle* particle,
	122	const AliMCParticle* mother,
	123	Int_t longest,
	124	Double_t vz,
	125	UShort_t nC,
	126	UShort_t nT,
2b556440	127	AliESDFMD& output,
2b556440	128	Double_t w) const
f68d9069	129	{
	130	// Store a particle.
	131	if (longest < 0) return;
	132
	133	AliTrackReference* ref = particle->GetTrackReference(longest);
	134	if (!ref) return;
	135
	136	// Get the detector coordinates
	137	UShort_t d, s, t;
	138	Char_t r;
	139	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
	140
	141	// Check if we have value already
	142	Double_t old = output.Multiplicity(d,r,s,t);
	143
	144	// If invalid, force it valid
	145	if (old == AliESDFMD::kInvalidMult) old = 0;
	146
	147	// Increment count
2b556440	148	output.SetMultiplicity(d,r,s,t,old+w);
f68d9069	149
	150	// Get track-reference stuff
	151	Double_t x = ref->X();
	152	Double_t y = ref->Y();
	153	Double_t z = ref->Z()-vz;
	154	Double_t rr = TMath::Sqrt(xx+yy);
	155	Double_t thetaR = TMath::ATan2(rr,z);
	156	Double_t phiR = TMath::ATan2(y,x);
	157	Double_t etaR = -TMath::Log(TMath::Tan(thetaR/2));
	158
	159	// Correct angle and convert to degrees
	160	if (thetaR < 0) thetaR += 2*TMath::Pi();
	161	thetaR *= 180. / TMath::Pi();
	162	if (phiR < 0) phiR += 2*TMath::Pi();
	163	phiR *= 180. / TMath::Pi();
	164
	165	// Fill histograms
	166	fNr->Fill(nC);
	167	fNt->Fill(nT);
	168
	169	const AliMCParticle* mp = (mother ? mother : particle);
	170	Double_t dEta = mp->Eta() - etaR;
	171	Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
	172	if (dPhi > 180) dPhi -= 360;
	173	if (dPhi < -180) dPhi += 360;
	174	fBinFlow->Fill(dEta, dPhi);
	175	fEtaBinFlow->Fill(etaR, dEta);
	176	fPhiBinFlow->Fill(phiR, dPhi);
	177
	178	if (fDebug)
	179	Info("StoreParticle", "Store @ FMD%d%c[%2d,%3d] nStrips=%3d, "
	180	"dEta=%7.4f, dPhi=%d",
	181	d, r, s, t, nT, dEta, int(dPhi+.5));
	182	}
	183
	184	//____________________________________________________________________
	185	Double_t
	186	AliFMDMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref,
	187	Double_t vz) const
	188	{
	189	// Get the incidient angle of the track reference.
	190	Double_t x = ref->X();
	191	Double_t y = ref->Y();
	192	Double_t z = ref->Z()-vz;
	193	Double_t rr = TMath::Sqrt(xx+yy);
	194	Double_t theta= TMath::ATan2(rr,z);
	195	Double_t ang = TMath::Abs(TMath::Pi()-theta);
	196	return ang;
	197	}
	198
	199	//____________________________________________________________________
	200	const AliMCParticle*
	201	AliFMDMCTrackDensity::GetMother(Int_t iTr,
	202	const AliMCEvent& event) const
	203	{
	204	//
	205	// Track down primary mother
	206	//
	207	Int_t i = iTr;
	208	do {
	209	const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
	210	if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
	211
	212	i = p->GetMother();
213	} while (i > 0);
214
215	return 0;
216	}
217
2b556440	218	// #define USE_FLOW_WEIGHTS 1
f68d9069	219	//____________________________________________________________________
	220	Bool_t
	221	AliFMDMCTrackDensity::Calculate(const AliESDFMD& input,
	222	const AliMCEvent& event,
	223	Double_t vz,
	224	AliESDFMD& output,
	225	TH2D* primary)
	226	{
	227	//
	228	// Filter the input kinematics and track references, using
	229	// some of the ESD information
	230	//
	231	// Parameters:
	232	// input Input ESD event
	233	// event Input MC event
	234	// vz Vertex position
	235	// output Output ESD-like object
	236	// primary Per-event histogram of primaries
	237	//
	238	// Return:
	239	// True on succes, false otherwise
	240	//
	241	output.Clear();
	242
	243	// Copy eta values to output
	244	for (UShort_t ed = 1; ed <= 3; ed++) {
	245	UShort_t nq = (ed == 1 ? 1 : 2);
	246	for (UShort_t eq = 0; eq < nq; eq++) {
	247	Char_t er = (eq == 0 ? 'I' : 'O');
	248	UShort_t ns = (eq == 0 ? 20 : 40);
	249	UShort_t nt = (eq == 0 ? 512 : 256);
	250	for (UShort_t es = 0; es < ns; es++)
	251	for (UShort_t et = 0; et < nt; et++)
	252	output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
	253	}
	254	}
2b556440	255
	256	#ifdef USE_FLOW_WEIGHTS
	257	AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
	258	(event.GenEventHeader());
	259	Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
	260	Double_t b = (hd ? hd->ImpactParameter() : -1 );
	261	#endif
	262
f68d9069	263	AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
	264	Int_t nTracks = stack->GetNtrack();//event.GetNumberOfTracks();
	265	Int_t nPrim = stack->GetNprimary();//event.GetNumberOfPrimary();
	266	for (Int_t iTr = 0; iTr < nTracks; iTr++) {
	267	AliMCParticle* particle =
	268	static_cast<AliMCParticle*>(event.GetTrack(iTr));
	269
	270	// Check the returned particle
	271	if (!particle) continue;
	272
	273	// Check if this charged and a primary
	274	Bool_t isCharged = particle->Charge() != 0;
	275	if (!isCharged) continue;
	276	Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
	277
	278	// Pseudo rapidity and azimuthal angle
	279	Double_t eta = particle->Eta();
	280	Double_t phi = particle->Phi();
	281
	282	// Fill 'dn/deta' histogram
	283	if (isPrimary && iTr < nPrim) {
	284	if (primary) primary->Fill(eta, phi);
	285	}
	286
	287	// Bail out if we're only processing primaries - perhaps we should
	288	// track back to the original primary?
	289	if (fUseOnlyPrimary && !isPrimary) continue;
	290
	291	Int_t nTrRef = particle->GetNumberOfTrackReferences();
	292	Int_t longest = -1;
	293	Double_t angle = 0;
	294	UShort_t oD = 0, oS = 1024, oT = 1024, ooT=1024;
	295	Char_t oR = '\0';
	296	UShort_t nC = 0;
	297	UShort_t nT = 0;
b4db5d3b	298	// Double_t oTheta= 0;
f68d9069	299	for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
	300	AliTrackReference* ref = particle->GetTrackReference(iTrRef);
	301
	302	// Check existence
	303	if (!ref) continue;
	304
	305	// Check that we hit an FMD element
	306	if (ref->DetectorId() != AliTrackReference::kFMD)
	307	continue;
	308
	309	// Get the detector coordinates
	310	UShort_t d, s, t;
	311	Char_t r;
	312	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
	313
	314	// Calculate length of last and second to last strips.
	315
	316	// IF base of cluster not set, set it here.
	317	if (ooT == 1024) ooT = t;
	318
	319	// Calculate distance of previous reference to base of cluster
	320	nT = TMath::Abs(t - ooT) + 1;
	321
	322	// Count number of track refs in this sector
	323	nC++;
	324
	325	Bool_t used = false;
	326	// If this is a new detector/ring, then reset the other one
	327	// Check if we have a valid old detectorm ring, and sector
	328	if (oD > 0 && oR != '\0' && oS != 1024) {
	329	// New detector, new ring, or new sector
	330	if (d != oD \|\| r != oR \|\| s != oS) {
	331	if (fDebug) Info("Process", "New because new sector");
	332	used = true;
	333	}
	334	}
	335	if (nT > fMaxConsequtiveStrips) {
	336	if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
	337	nT, t, oT, ooT);
	338	used = true;
	339	}
	340	if (used) {
	341	if (fDebug)
	342	Info("Process", "I=%3d L=%3d D=%d (was %d), R=%c (was %c), "
	343	"S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
	344	iTr, longest, d, oD, r, oR, s, oS, t, oT, nT,
	345	fMaxConsequtiveStrips);
	346	Int_t nnT = TMath::Abs(oT - ooT) + 1;
	347	const AliMCParticle* mother = GetMother(iTr, event);
2b556440	348
	349	#ifdef USE_FLOW_WEIGHTS
	350	phi = (mother ? mother->Phi() : particle->Phi());
	351	eta = (mother ? mother->Eta() : particle->Eta());
	352	Double_t pt = (mother ? mother->Pt() : particle->Pt());
	353	Int_t id = (mother ? mother->PdgCode() : 2212);
	354	Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
	355	#else
	356	Double_t weight = 1; // //
	357	#endif
	358	StoreParticle(particle, mother, longest, vz, nC, nnT, output, weight);
f68d9069	359	longest = -1;
	360	angle = 0;
	361	nC = 1; // Reset track-ref counter - we have this->1
	362	nT = 0; // Reset cluster size
	363	oD = 0; // Reset detector
	364	oR = '\0'; // Reset ring
	365	oS = 1024; // Reset sector
	366	oT = 1024; // Reset old strip
	367	ooT = t; // Reset base
	368	}
	369	else if (fDebug) {
	370	if (ooT == t)
	371	Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
	372	d, r, s, t);
	373	else
	374	Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
	375	"length=%3d (now in %3d, previous %3d)",
	376	d, r, s, ooT, nT, t, oT);
	377	}
	378
	379	oD = d;
	380	oR = r;
	381	oS = s;
	382	oT = t;
	383	nT = TMath::Abs(t-ooT)+1;
	384
	385	// The longest passage is determined through the angle
	386	Double_t ang = GetTrackRefTheta(ref, vz);
	387	if (ang > angle) {
	388	longest = iTrRef;
	389	angle = ang;
	390	}
b4db5d3b	391	// oTheta = ang;
f68d9069	392	} // Loop over track references
	393	if (longest < 0) continue; // Nothing found
	394
	395	// Get the reference corresponding to the longest path through the detector
	396	if (fDebug)
	397	Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)",
	398	iTr, longest, nT, fMaxConsequtiveStrips);
	399	const AliMCParticle* mother = GetMother(iTr, event);
2b556440	400
	401	#ifdef USE_FLOW_WEIGHTS
	402	phi = (mother ? mother->Phi() : particle->Phi());
	403	eta = (mother ? mother->Eta() : particle->Eta());
	404	Double_t pt = (mother ? mother->Pt() : particle->Pt());
	405	Int_t id = (mother ? mother->PdgCode() : 2212);
	406	Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
	407	#else
	408	Double_t weight = 1; // //
	409	#endif
	410	StoreParticle(particle, mother, longest, vz, nC, nT, output, weight);
f68d9069	411	} // Loop over tracks
	412	return kTRUE;
	413	}
2b556440	414
	415	//____________________________________________________________________
	416	Double_t
	417	AliFMDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b,
	418	Double_t phi, Double_t rp, Int_t id) const
	419	{
	420	static TF1 gaus = TF1("gaus", "gaus", -6, 6);
	421	gaus.SetParameters(0.1, 0., 9);
	422	// gaus.SetParameters(0.1, 0., 3);
	423	// gaus.SetParameters(0.1, 0., 15);
	424
	425	const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350,
	426	0.45, 0.55, 0.650,
	427	0.75, 0.85, 0.950,
	428	1.10, 1.30, 1.500,
	429	1.70, 1.90, 2.250,
	430	2.75, 3.25, 3.750,
	431	4.50};
	432	const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400,
	433	0.059911,0.073516,
	434	0.089756,0.105486,
	435	0.117391,0.128199,
	436	0.138013,0.158271,
	437	0.177726,0.196383,
	438	0.208277,0.216648,
	439	0.242954,0.249961,
	440	0.240131,0.269006,
	441	0.207796};
	442	const Int_t nPointsCumulant2nd4050ALICE =
	443	sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);
	444	static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
	445	#if 0
	446	const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35,
	447	0.45,0.550,0.65,
	448	0.75,0.850,0.95,
	449	1.10,1.300,1.50,
	450	1.70,1.900,2.25,
	451	2.75,3.250,3.75,
	452	4.50,5.500,7.00,
	453	9.000000};
	454	const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,
	455	0.048566,0.061083,
	456	0.070910,0.078831,
	457	0.091396,0.102026,
	458	0.109691,0.124449,
	459	0.139819,0.155561,
	460	0.165701,0.173678,
	461	0.191149,0.202015,
	462	0.204540,0.212560,
	463	0.195885,0.000000,
	464	0.000000,0.000000};
	465	#endif
	466	const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350,
	467	0.45,0.55,0.650,
	468	0.75,0.85,0.950,
	469	1.10,1.30,1.500,
	470	1.70,1.90,2.250,
	471	2.75,3.25,3.750,
	472	4.50};
	473	const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,
	474	0.049824,0.066662,
	475	0.075856,0.081583,
	476	0.099778,0.104674,
	477	0.118545,0.131874,
478	0.152959,0.155348,
479	0.169751,0.179052,
480	0.178532,0.198851,
481	0.185737,0.239901,
482	0.186098};
483	const Int_t nPointsCumulant4th4050ALICE =
484	sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);
485	static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE,
486	xCumulant4th4050ALICE,
487	yCumulant4th4050ALICE);
488
489	const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75,
490	4.225,
491	5.965,
492	7.765,
493	9.215,
494	10.46,
495	11.565,
496	12.575};
497	const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,
498	0.055818,0.073137,
499	0.083898,0.086690,
500	0.082040,0.077777};
501	const Int_t nPointsCumulant4thTPCrefMultTPConlyAll =
502	sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
503	TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll,
504	xCumulant4thTPCrefMultTPConlyAll,
505	yCumulant4thTPCrefMultTPConlyAll);
506
507
508	Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 +
509	alicePointsPt4.Eval(pt) * 0.5)
510	* (aliceCent.Eval(b) / aliceCent.Eval(10.46))
511	* 2. * TMath::Cos(2. * (phi - rp)));
512	if (TMath::Abs(id) == 211) weight *= 1.3; //pion flow
513	else if (TMath::Abs(id) == 2212) weight *= 1.0; //proton flow
514	else weight *= 0.7;
515
516	return weight;
517	}
f68d9069	518	//____________________________________________________________________
	519	void
	520	AliFMDMCTrackDensity::Print(Option_t* /option/) const
	521	{
	522	char ind[gROOT->GetDirLevel()+1];
	523	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	524	ind[gROOT->GetDirLevel()] = '\0';
	525	std::cout << ind << ClassName() << ": " << GetName() << '\n'
	526	<< std::boolalpha
	527	<< ind << " Only primary tracks: " << fUseOnlyPrimary << '\n'
	528	<< ind << " Max cluster size: " << fMaxConsequtiveStrips
	529	<< std::noboolalpha << std::endl;
	530
	531	}
	532
	533	//____________________________________________________________________
	534	//
	535	// EOF
	536	//