PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChainESD.cxx

   1 /////////////////////////////////////////////////////////////////////////////////////
   2 //                                                                                 //
   3 // AliFemtoEventReaderKinematicsChainESD - the reader class for the Alice ESD and     //
   4 // the model Kinematics information tailored for the Task framework and the        //
   5 // Reads in AliESDfriend to create shared hit/quality information                  //
   6 // Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
   7 //          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
   8 //                                                                                 //
   9 /////////////////////////////////////////////////////////////////////////////////////
  10
  11 #include "AliFemtoEventReaderKinematicsChainESD.h"
  12
  13 #include "TFile.h"
  14 #include "TTree.h"
  15 #include "TList.h"
  16
  17 #include "AliFmPhysicalHelixD.h"
  18 #include "AliFmThreeVectorF.h"
  19
  20 #include "SystemOfUnits.h"
  21
  22 #include "AliFemtoEvent.h"
  23
  24 #include "TParticle.h"
  25 #include "AliStack.h"
  26 #include "AliESDEvent.h"
  27 #include "AliCentrality.h"
  28 #include "TParticlePDG.h"
  29 #include "AliFemtoModelHiddenInfo.h"
  30 #include "AliFemtoModelGlobalHiddenInfo.h"
  31 #include "AliGenHijingEventHeader.h"
  32 #include "AliGenCocktailEventHeader.h"
  33
  34 #include "AliVertexerTracks.h"
  35
  36 ClassImp(AliFemtoEventReaderKinematicsChainESD)
  37
  38 #if !(ST_NO_NAMESPACES)
  39   using namespace units;
  40 #endif
  41
  42 using namespace std;
  43 //____________________________
  44 AliFemtoEventReaderKinematicsChainESD::AliFemtoEventReaderKinematicsChainESD():
  45   fFileName(" "),
  46   fConstrained(true),
  47   fNumberofEvent(0),
  48   fCurEvent(0),
  49   fCurFile(0),
  50   fStack(0x0),
  51   fEvent(0x0),
  52   fGenHeader(0x0),
  53   fEstEventMult(kGlobalCount),
  54   fRotateToEventPlane(0)
  55 {
  56   //constructor with 0 parameters , look at default settings
  57 }
  58
  59 //__________________
  60 AliFemtoEventReaderKinematicsChainESD::AliFemtoEventReaderKinematicsChainESD(const AliFemtoEventReaderKinematicsChainESD& aReader):
  61   AliFemtoEventReader(aReader),
  62   fFileName(" "),
  63   fConstrained(true),
  64   fNumberofEvent(0),
  65   fCurEvent(0),
  66   fCurFile(0),
  67   fStack(0x0),
  68   fEvent(0x0),
  69   fGenHeader(0x0),
  70   fEstEventMult(kGlobalCount),
  71   fRotateToEventPlane(0)
  72 {
  73   // Copy constructor
  74   fConstrained = aReader.fConstrained;
  75   fNumberofEvent = aReader.fNumberofEvent;
  76   fCurEvent = aReader.fCurEvent;
  77   fCurFile = aReader.fCurFile;
  78   fStack = aReader.fStack;
  79   fEvent = aReader.fEvent;
  80   fEstEventMult = aReader.fEstEventMult;
  81   fRotateToEventPlane = aReader.fRotateToEventPlane;
  82 }
  83 //__________________
  84 AliFemtoEventReaderKinematicsChainESD::~AliFemtoEventReaderKinematicsChainESD()
  85 {
  86   //Destructor
  87   //delete fEvent;
  88 }
  89
  90 //__________________
  91 AliFemtoEventReaderKinematicsChainESD& AliFemtoEventReaderKinematicsChainESD::operator=(const AliFemtoEventReaderKinematicsChainESD& aReader)
  92 {
  93   // Assignment operator
  94   if (this == &aReader)
  95     return *this;
  96
  97   fConstrained = aReader.fConstrained;
  98   fNumberofEvent = aReader.fNumberofEvent;
  99   fCurEvent = aReader.fCurEvent;
 100   fCurFile = aReader.fCurFile;
 101   fStack = aReader.fStack;
 102   fEvent = aReader.fEvent;
 103   fGenHeader = aReader.fGenHeader;
 104   fEstEventMult = aReader.fEstEventMult;
 105   fRotateToEventPlane = aReader.fRotateToEventPlane;
 106   return *this;
 107 }
 108 //__________________
 109 // Simple report
 110 AliFemtoString AliFemtoEventReaderKinematicsChainESD::Report()
 111 {
 112   AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChainESD\n";
 113   return temp;
 114 }
 115
 116 //__________________
 117 void AliFemtoEventReaderKinematicsChainESD::SetConstrained(bool constrained)
 118 {
 119   // Select whether to read constrained or not constrained momentum
 120   fConstrained=constrained;
 121 }
 122 //__________________
 123 bool AliFemtoEventReaderKinematicsChainESD::GetConstrained() const
 124 {
 125   // Check whether we read constrained or not constrained momentum
 126   return fConstrained;
 127 }
 128 //__________________
 129 AliFemtoEvent* AliFemtoEventReaderKinematicsChainESD::ReturnHbtEvent()
 130 {
 131   // Get the event, read all the relevant information from the stack
 132   // and fill the AliFemtoEvent class
 133   // Returns a valid AliFemtoEvent
 134   AliFemtoEvent *hbtEvent = 0;
 135   string tFriendFileName;
 136
 137   cout << "AliFemtoEventReaderKinematicsChainESD::Starting to read event: "<<fCurEvent<<endl;
 138
 139   hbtEvent = new AliFemtoEvent;
 140   //setting basic things
 141   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 142   hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
 143   hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
 144   hbtEvent->SetZDCN1Energy(0);
 145   hbtEvent->SetZDCP1Energy(0);
 146   hbtEvent->SetZDCN2Energy(0);
 147   hbtEvent->SetZDCP2Energy(0);
 148   hbtEvent->SetZDCEMEnergy(0);
 149   hbtEvent->SetZDCParticipants(0);
 150   hbtEvent->SetTriggerMask(0);
 151   hbtEvent->SetTriggerCluster(0);
 152
 153   //Vertex
 154   double fV1[3] = {0.0,0.0,0.0};
 155   double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 156
 157
 158   AliFmThreeVectorF vertex(0,0,0);
 159
 160
 161   hbtEvent->SetPrimVertPos(vertex);
 162   hbtEvent->SetPrimVertCov(fVCov);
 163
 164   Double_t tReactionPlane = 0;
 165
 166   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 167   if (!hdh) {
 168     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 169     if (cdh) {
 170       TList *tGenHeaders = cdh->GetHeaders();
 171       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 172         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 173         if (hdh) break;
 174       }
 175     }
 176   }
 177
 178   if (hdh)
 179     {
 180       tReactionPlane = hdh->ReactionPlaneAngle();
 181       cout << "Got reaction plane " << tReactionPlane << endl;
 182     }
 183
 184   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 185
 186   //starting to reading tracks
 187   int nofTracks=0;  //number of all tracks in MC event
 188   nofTracks=fStack->GetNtrack();
 189   int realnofTracks=0;//number of track which we use in analysis
 190
 191
 192   int tNormMult = 0;
 193   int tV0direction = 0;
 194   for (int i=0;i<nofTracks;i++)
 195     {
 196       //take only primaries
 197       if(!fStack->IsPhysicalPrimary(i)) {continue;}
 198
 199       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 200
 201           //getting next track
 202       TParticle *kinetrack= fStack->Particle(i);
 203
 204       //setting multiplicity
 205         realnofTracks++;//real number of tracks (only primary particles)
 206
 207       //setting normalized multiplicity
 208
 209         if(kinetrack->GetPDG()->Charge()/3!=0)
 210           if (kinetrack->Pt() > 0.15 && kinetrack->Pt() < 20)
 211             if (kinetrack->Eta() < 0.8)
 212               tNormMult++;
 213
 214         //counting particles that go into direction of VZERO detector
 215         if(kinetrack->Eta() > 2.8 && kinetrack->Eta() < 5.1) //VZERO-A
 216           tV0direction++;
 217         if(kinetrack->Eta() > -3.7 && kinetrack->Eta() < -1.7)//VZERO-C
 218           tV0direction++;
 219
 220           //charge
 221       trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
 222
 223
 224       //in aliroot we have AliPID
 225       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 226       //we use only 5 first
 227       double kinepid[5];
 228       for(int pid_iter=0;pid_iter<5;pid_iter++)
 229           kinepid[pid_iter]=0;
 230
 231       int pdgcode = kinetrack->GetPdgCode();
 232       //proton
 233       if(pdgcode==2212 || pdgcode==-2212)
 234         kinepid[4]=1000;
 235       //kaon
 236       else if(pdgcode==321 || pdgcode==-321 )
 237         kinepid[3]=1000;
 238       //pion
 239       else if( pdgcode==211 || pdgcode==-211)
 240         kinepid[2]=1000;
 241       //electron
 242       else if(pdgcode==11 || pdgcode==-11)
 243         kinepid[0]=1000;
 244       //muon
 245       else if(pdgcode==13 || pdgcode==-13)
 246         kinepid[1]=1000;
 247       else
 248         continue;
 249       trackCopy->SetPidProbElectron(kinepid[0]);
 250       trackCopy->SetPidProbMuon(kinepid[1]);
 251       trackCopy->SetPidProbPion(kinepid[2]);
 252       trackCopy->SetPidProbKaon(kinepid[3]);
 253       trackCopy->SetPidProbProton(kinepid[4]);
 254
 255
 256         //Momentum
 257       double pxyz[3];
 258       // double rxyz[3];
 259
 260         pxyz[0]=kinetrack->Px();
 261         pxyz[1]=kinetrack->Py();
 262         pxyz[2]=kinetrack->Pz();
 263
 264         // rxyz[0]=kinetrack->Vx();
 265         // rxyz[1]=kinetrack->Vy();
 266         // rxyz[2]=kinetrack->Vz();
 267
 268         if (fRotateToEventPlane) {
 269           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 270           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 271
 272           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 273           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 274         }
 275
 276         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 277         if (v.Mag() < 0.0001) {
 278           //cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 279           delete trackCopy;
 280           continue;
 281         }
 282
 283         trackCopy->SetP(v);//setting momentum
 284         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 285         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 286         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 287
 288         //label
 289         trackCopy->SetLabel(i);
 290
 291
 292         hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 293         //cout<<"Track added: "<<i<<endl;
 294
 295     }
 296
 297   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 298
 299   if (fEstEventMult == kGlobalCount)
 300     hbtEvent->SetNormalizedMult(tNormMult);
 301   else if(fEstEventMult == kVZERO)
 302     hbtEvent->SetNormalizedMult(tV0direction);
 303   else if (fEstEventMult == kReferenceITSTPC)
 304     {
 305       hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.2));
 306     }
 307   else if (fEstEventMult == kCentrality) {
 308     // centrality between 0 (central) and 1 (very peripheral)
 309     AliCentrality *cent = fEvent->GetCentrality();
 310     if (cent) {
 311       if (cent->GetCentralityPercentile("V0M") < 0.00001)
 312         hbtEvent->SetNormalizedMult(-1);
 313       else
 314         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
 315       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 316                              10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
 317     }
 318   }
 319
 320   fCurEvent++;
 321
 322
 323   //V0 analysis code - no V0 finder for Kinematics, we can only check if it is primary and if it has at least 2 daughters.
 324
 325   for (int i=0;i<nofTracks;i++)
 326     {
 327       //do not take primaries
 328       if(!fStack->IsPhysicalPrimary(i)) {continue;}
 329       //getting next track
 330       TParticle *kinetrack= fStack->Particle(i);
 331       if (!kinetrack) continue;
 332
 333       if(kinetrack->GetPDG()->Charge()!=0) continue; //charge - neutral
 334       //if(kinetrack->GetPDG()->Stable()==1) continue; //particle is not stable
 335       if(kinetrack->GetDaughter(0)<1) continue; //has 1'st daughter
 336       if(kinetrack->GetDaughter(1)<1) continue;  //has 2'nd daughter
 337
 338       //we want one positive, one negative particle. Or two neutral.
 339       // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()>=0)
 340       //        if((fStack->Particle(kinetrack->GetDaughter(1)))->GetPDG()->Charge()>0)
 341       //          continue;
 342       // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<=0)
 343       //        if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<0)
 344       //          continue;
 345
 346       if(kinetrack->Pt()<0.00001)
 347         continue;
 348
 349       AliFemtoV0* trackCopyV0 = new AliFemtoV0();
 350       CopyAODtoFemtoV0(kinetrack, trackCopyV0);
 351       hbtEvent->V0Collection()->push_back(trackCopyV0);
 352     //cout<<"Pushback v0 to v0collection"<<endl;
 353     }
 354
 355
 356   cout<<"Number of tracks: "<<realnofTracks<<endl;
 357
 358   return hbtEvent;
 359 }
 360
 361 //___________________
 362 void AliFemtoEventReaderKinematicsChainESD::SetStackSource(AliStack *aStack)
 363 {
 364   // The chain loads the stack for us
 365   // You must provide the address where it can be found
 366   fStack = aStack;
 367 }
 368 void AliFemtoEventReaderKinematicsChainESD::SetESDSource(AliESDEvent *aESD)
 369 {
 370   // The chain loads the ESD for us
 371   // You must provide the address where it can be found
 372   fEvent = aESD;
 373 }
 374 //___________________
 375 void AliFemtoEventReaderKinematicsChainESD::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 376 {
 377   // The chain loads the generator event header for us
 378   // You must provide the address where it can be found
 379   fGenHeader = aGenHeader;
 380 }
 381
 382 //__________________
 383 void AliFemtoEventReaderKinematicsChainESD::SetRotateToEventPlane(short dorotate)
 384 {
 385   fRotateToEventPlane=dorotate;
 386 }
 387
 388 void AliFemtoEventReaderKinematicsChainESD::SetUseMultiplicity(EstEventMult aType)
 389 {
 390   fEstEventMult = aType;
 391 }
 392
 393 Float_t AliFemtoEventReaderKinematicsChainESD::GetSigmaToVertex(double *impact, double *covar)
 394 {
 395   // Calculates the number of sigma to the vertex.
 396
 397   Float_t b[2];
 398   Float_t bRes[2];
 399   Float_t bCov[3];
 400
 401   b[0] = impact[0];
 402   b[1] = impact[1];
 403   bCov[0] = covar[0];
 404   bCov[1] = covar[1];
 405   bCov[2] = covar[2];
 406
 407   bRes[0] = TMath::Sqrt(bCov[0]);
 408   bRes[1] = TMath::Sqrt(bCov[2]);
 409
 410   // -----------------------------------
 411   // How to get to a n-sigma cut?
 412   //
 413   // The accumulated statistics from 0 to d is
 414   //
 415   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 416   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 417   //
 418   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 419   // Can this be expressed in a different way?
 420
 421   if (bRes[0] == 0 || bRes[1] ==0)
 422     return -1;
 423
 424   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 425
 426   // stupid rounding problem screws up everything:
 427   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 428   if (TMath::Exp(-d * d / 2) < 1e-10)
 429     return 1000;
 430
 431   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 432   return d;
 433 }
 434
 435
 436
 437  void AliFemtoEventReaderKinematicsChainESD::CopyAODtoFemtoV0(TParticle *tv0, AliFemtoV0 *tFemtoV0 )
 438 {
 439   tFemtoV0->SetEtaV0(tv0->Eta());
 440   tFemtoV0->SetEtaV0(tv0->Phi());
 441   tFemtoV0->SetptV0(tv0->Pt());
 442   tFemtoV0->SetptotV0(tv0->P());
 443
 444   tFemtoV0->SetmomV0X(tv0->Px());
 445   tFemtoV0->SetmomV0Y(tv0->Py());
 446   tFemtoV0->SetmomV0Z(tv0->Pz());
 447   AliFemtoThreeVector momv0(tv0->Px(),tv0->Py(),tv0->Pz());
 448   tFemtoV0->SetmomV0(momv0);
 449
 450
 451   TParticle *trackpos;
 452   TParticle *trackneg;
 453
 454   //daughters
 455   if(fStack->Particle(tv0->GetDaughter(0))->GetPDG()->Charge()>=0) //first positive, second negative
 456     {
 457       trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
 458       trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
 459       tFemtoV0->SetidPos(tv0->GetDaughter(0));
 460       tFemtoV0->SetidNeg(tv0->GetDaughter(1));
 461     }
 462   else //first negative, second positive
 463     {
 464       trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
 465       trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
 466       tFemtoV0->SetidPos(tv0->GetDaughter(1));
 467       tFemtoV0->SetidNeg(tv0->GetDaughter(0));
 468     }
 469
 470   tFemtoV0->SetEtaPos(trackpos->Eta());
 471   tFemtoV0->SetEtaNeg(trackneg->Eta());
 472
 473   tFemtoV0->SetptPos(trackpos->Pt());
 474   tFemtoV0->SetptNeg(trackneg->Pt());
 475
 476   tFemtoV0->SetptotPos(trackpos->P());
 477   tFemtoV0->SetptotNeg(trackneg->P());
 478
 479   tFemtoV0->SetmomPosX(trackpos->Px());
 480   tFemtoV0->SetmomPosY(trackpos->Py());
 481   tFemtoV0->SetmomPosZ(trackpos->Pz());
 482   AliFemtoThreeVector mompos(trackpos->Px(),trackpos->Py(),trackpos->Pz());
 483   tFemtoV0->SetmomPos(mompos);
 484
 485   tFemtoV0->SetmomNegX(trackneg->Px());
 486   tFemtoV0->SetmomNegY(trackneg->Py());
 487   tFemtoV0->SetmomNegZ(trackneg->Pz());
 488   AliFemtoThreeVector momneg(trackneg->Px(),trackneg->Py(),trackneg->Pz());
 489   tFemtoV0->SetmomNeg(momneg);
 490
 491
 492   tFemtoV0->SetmassLambda(tv0->GetMass());
 493   tFemtoV0->SetmassAntiLambda(tv0->GetMass());
 494   tFemtoV0->SetmassK0Short(tv0->GetMass());
 495
 496   tFemtoV0->SetYV0(tv0->Y());
 497
 498   tFemtoV0->SetdecayVertexV0X(trackpos->Vx()); //vertex of the decay is set as the vertex of creation of daughters
 499   tFemtoV0->SetdecayVertexV0Y(trackpos->Vy());
 500   tFemtoV0->SetdecayVertexV0Z(trackpos->Vz());
 501   AliFemtoThreeVector decayvertex(trackpos->Vx(),trackpos->Vy(),trackpos->Vz());
 502   tFemtoV0->SetdecayVertexV0(decayvertex);
 503
 504   tFemtoV0->SetdcaV0Daughters(0);
 505   tFemtoV0->SetCosPointingAngle(1);
 506
 507
 508   tFemtoV0->SetStatusPos(1);
 509   tFemtoV0->SetStatusNeg(1);
 510
 511
 512   if(trackpos->GetPdgCode()==2212) //proton
 513     {
 514       tFemtoV0->SetPosNSigmaTPCK(1000);
 515       tFemtoV0->SetPosNSigmaTPCPi(1000);
 516       tFemtoV0->SetPosNSigmaTPCP(0);
 517     }
 518   if(trackneg->GetPdgCode()==-2212) //antiproton
 519     {
 520       tFemtoV0->SetNegNSigmaTPCK(1000);
 521       tFemtoV0->SetNegNSigmaTPCPi(1000);
 522       tFemtoV0->SetNegNSigmaTPCP(0);
 523     }
 524   if(trackpos->GetPdgCode()==211) //pion plus
 525     {
 526       tFemtoV0->SetPosNSigmaTPCK(1000);
 527       tFemtoV0->SetPosNSigmaTPCPi(0);
 528       tFemtoV0->SetPosNSigmaTPCP(1000);
 529     }
 530   if(trackneg->GetPdgCode()==-211) //pion minus
 531     {
 532       tFemtoV0->SetNegNSigmaTPCK(1000);
 533       tFemtoV0->SetNegNSigmaTPCPi(0);
 534       tFemtoV0->SetNegNSigmaTPCP(1000);
 535     }
 536   if(trackpos->GetPdgCode()==321) //K+
 537     {
 538       tFemtoV0->SetPosNSigmaTPCK(0);
 539       tFemtoV0->SetPosNSigmaTPCPi(1000);
 540       tFemtoV0->SetPosNSigmaTPCP(1000);
 541     }
 542   if(trackneg->GetPdgCode()==-321) //K-
 543     {
 544       tFemtoV0->SetNegNSigmaTPCK(0);
 545       tFemtoV0->SetNegNSigmaTPCPi(1000);
 546       tFemtoV0->SetNegNSigmaTPCP(1000);
 547     }
 548
 549
 550 }