PWGCF/FEMTOSCOPY/AliFemto/AliFemtoSpherocityEventCut.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoSpherocityEventCut - the basic cut for events.                     //
   4 // Only cuts on event multiplicity, z-vertex position and                     //
   5 // transverse spherocity are accepted.                                        //
   6 //                                                                            //
   7 ////////////////////////////////////////////////////////////////////////////////
   8
   9 #include "AliFemtoSpherocityEventCut.h"
  10 //#include <cstdio>
  11
  12 #ifdef __ROOT__
  13 ClassImp(AliFemtoSpherocityEventCut)
  14 #endif
  15
  16 AliFemtoSpherocityEventCut::AliFemtoSpherocityEventCut() :
  17   AliFemtoEventCut(),
  18   fEventMult(),
  19   fVertZPos(),
  20   fAcceptBadVertex(false),
  21   fNEventsPassed(0),
  22   fNEventsFailed(0),
  23   fAcceptOnlyPhysics(0),
  24   fSoCutMin(0.0),
  25   fSoCutMax(1.0),
  26   fSelectTrigger(0)
  27 {
  28   // Default constructor
  29   fEventMult[0] = 0;
  30   fEventMult[1] = 100000;
  31   fVertZPos[0] = -100.0;
  32   fVertZPos[1] = 100.0;
  33   fPsiEP[0] = -1000.0;
  34   fPsiEP[1] = 1000.0;
  35 }
  36 //------------------------------
  37 AliFemtoSpherocityEventCut::~AliFemtoSpherocityEventCut(){
  38   // Default destructor
  39 }
  40 //------------------------------
  41 bool AliFemtoSpherocityEventCut::Pass(const AliFemtoEvent* event){
  42
  43   // Pass events if they fall within the multiplicity and z-vertex
  44   // position range. Fail otherwise
  45   //  int mult =  event->NumberOfTracks();
  46   int mult = (int) event->UncorrectedNumberOfPrimaries();
  47   double vertexZPos = event->PrimVertPos().z();
  48   double spherocity=-10;
  49   int MULT=0;
  50
  51
  52    AliFemtoTrackCollection * tracks = event->TrackCollection();
  53   for (AliFemtoTrackIterator iter=tracks->begin();iter!=tracks->end();iter++){
  54
  55     Double_t NewPt =  (*iter)->Pt();
  56     Double_t NewEta = (*iter)->P().PseudoRapidity();
  57     if(TMath::Abs(NewEta)>0.8 || NewPt<0.5){continue;}
  58
  59      MULT++;
  60
  61   }
  62     //if(SumPt==0){return kFALSE;}
  63       if(MULT<3){return kFALSE;}
  64
  65         Double_t *pxA=new Double_t[MULT]();
  66         Double_t *pyA=new Double_t[MULT]();
  67
  68         Double_t sumapt=0;
  69         Int_t counter=0;
  70
  71    AliFemtoTrackCollection * tracks2 = event->TrackCollection();
  72   for (AliFemtoTrackIterator iter2=tracks2->begin();iter2!=tracks2->end();iter2++){
  73
  74     Double_t NewPt2 =  (*iter2)->Pt();
  75     Double_t NewPhi2 = (*iter2)->P().Phi();
  76     Double_t NewEta2 = (*iter2)->P().PseudoRapidity();
  77     if(TMath::Abs(NewEta2)>0.8 || NewPt2<0.5){continue;}
  78
  79     Double_t Px;
  80     Double_t Py;
  81
  82     Px= NewPt2 * TMath::Cos(NewPhi2);
  83     Py= NewPt2 * TMath::Sin(NewPhi2);
  84
  85
  86     pxA[counter]=Px;
  87     pyA[counter]=Py;
  88     sumapt+=NewPt2;
  89     counter++;
  90
  91   }
  92
  93   Double_t pFull = 0;
  94   Double_t Spherocity = 2;
  95   //Getting thrust
  96   for(Int_t i = 0; i < 360; ++i){
  97     Double_t numerador = 0;
  98     Double_t phiparam  = 0;
  99     Double_t nx = 0;
 100     Double_t ny = 0;
 101     phiparam=((TMath::Pi()) * i) / 180;   // parametrization of the angle
 102     nx = TMath::Cos(phiparam);            // x component of an unitary vector n
 103     ny = TMath::Sin(phiparam);            // y component of an unitary vector n
 104     for(Int_t i1 = 0; i1 < MULT; ++i1){
 105       numerador += TMath::Abs(ny * pxA[i1] - nx * pyA[i1]);//product between momentum proyection in XY plane and the unitari vector.
 106     }
 107     pFull=TMath::Power( (numerador / sumapt),2 );
 108     if(pFull < Spherocity)//maximization of pFull
 109       {Spherocity = pFull;}
 110   }
 111
 112
 113     spherocity=((Spherocity)*TMath::Pi()*TMath::Pi())/4.0;
 114
 115   if(pxA){// clean up array memory used for TMath::Sort
 116     delete[] pxA;
 117     pxA=0;
 118   }
 119   if(pyA){// clean up array memory used for TMath::Sort
 120     delete[] pyA;
 121     pyA=0;
 122   }
 123
 124   if(spherocity>fSoCutMax || spherocity<fSoCutMin){
 125           //cout<<" Event kicked out !"<<"SoCutMax= "<<fSoCutMax<<"  SoCutMin= "<<fSoCutMin<<endl;
 126           return kFALSE;}
 127
 128   double epvzero = event->ReactionPlaneAngle();
 129
 130   // cout << "AliFemtoSpherocityEventCut:: epvzero:       " << fPsiEP[0] << " < " << epvzero << " < " << fPsiEP[1] << endl;
 131 //   cout << "AliFemtoSpherocityEventCut:: mult:       " << fEventMult[0] << " < " << mult << " < " << fEventMult[1] << endl;
 132 //   cout << "AliFemtoSpherocityEventCut:: VertexZPos: " << fVertZPos[0] << " < " << vertexZPos << " < " << fVertZPos[1] << endl;
 133 //   cout << "AliFemtoSpherocityEventCut:: VertexZErr: " << event->PrimVertCov()[4] << endl;
 134
 135   // cout << "AliFemtoSpherocityEventCut:: MagneticField: " << event->MagneticField() << endl;
 136   // cout << "AliFemtoSpherocityEventCut:: IsCollisionCandidate: " << event->IsCollisionCandidate() << endl;
 137   // cout << "AliFemtoSpherocityEventCut:: TriggerCluster: " << event->TriggerCluster() << endl;
 138   // cout << "AliFemtoSpherocityEventCut:: fSelectTrigger: " << fSelectTrigger << endl;
 139   // cout << "AliFemtoSpherocityEventCut:: " << endl;
 140   bool goodEvent =
 141     ((mult >= fEventMult[0]) &&
 142      (mult <= fEventMult[1]) &&
 143      (vertexZPos > fVertZPos[0]) &&
 144      (vertexZPos < fVertZPos[1]) &&
 145      (epvzero > fPsiEP[0]) &&
 146      (epvzero < fPsiEP[1]) &&
 147      ((!fAcceptBadVertex) || (event->ZDCParticipants() > 1.0)) &&
 148       ((!fSelectTrigger) || (event->TriggerCluster() == fSelectTrigger))
 149 );
 150
 151   // cout << "AliFemtoSpherocityEventCut:: goodEvent" <<goodEvent << endl;
 152
 153   goodEvent ? fNEventsPassed++ : fNEventsFailed++ ;
 154   // cout << "AliFemtoSpherocityEventCut:: return : " << goodEvent << endl;
 155 //     (fAcceptBadVertex || (event->PrimVertCov()[4] > -1000.0)) &&
 156
 157   return (goodEvent);
 158 }
 159 //------------------------------
 160 AliFemtoString AliFemtoSpherocityEventCut::Report(){
 161   // Prepare report
 162   string stemp;
 163   char ctemp[100];
 164   snprintf(ctemp , 100, "\nMultiplicity:\t %d-%d",fEventMult[0],fEventMult[1]);
 165   stemp = ctemp;
 166   snprintf(ctemp , 100, "\nVertex Z-position:\t %E-%E",fVertZPos[0],fVertZPos[1]);
 167   stemp += ctemp;
 168   snprintf(ctemp , 100, "\nNumber of events which passed:\t%ld  Number which failed:\t%ld",fNEventsPassed,fNEventsFailed);
 169   stemp += ctemp;
 170   AliFemtoString returnThis = stemp;
 171   return returnThis;
 172 }
 173 void AliFemtoSpherocityEventCut::SetAcceptBadVertex(bool b)
 174 {
 175   fAcceptBadVertex = b;
 176 }
 177 bool AliFemtoSpherocityEventCut::GetAcceptBadVertex()
 178 {
 179   return fAcceptBadVertex;
 180 }