PWG2/FLOW/AliFlowWeighter.cxx

   1 //////////////////////////////////////////////////////////////////////
   2 //
   3 // $Id$
   4 //
   5 // Author: Emanuele Simili
   6 //
   7 //////////////////////////////////////////////////////////////////////
   8 //_____________________________________________________________
   9 //
  10 // Description:
  11 //         the AliFlowWeighter class generates the phi-weights which
  12 // are later used to correct for azimuthal anisotropy in the reconstruction
  13 // efficiency of the ALICE TPC. It also fills an histogram of normalised
  14 // particle abundancies, which can be used as bayesian weights for particle Id.
  15 //
  16 // - The method Init() generates the histograms and opens a new fPhiWgt file.
  17 // - The method WeightEvent(AliFlowEvent*) fills phi and PId histograms.
  18 //   It must be inserted in a loop over the event sample.
  19 // - The method Finish() calculates the weights, saves the histograms and
  20 //   closes the file. The AliFlowSelection object is saved as well.
  21 //
  22 //////////////////////////////////////////////////////////////////////
  23
  24 #ifndef ALIFLOWWEIGHTER_CXX
  25 #define ALIFLOWWEIGHTER_CXX
  26
  27 // ROOT things
  28 #include <TROOT.h>
  29 #include <TMath.h>
  30 #include <TFile.h>
  31 #include <TTree.h>
  32 #include <TString.h>
  33 #include <TObject.h>
  34 #include <TObjArray.h>
  35 #include <TOrdCollection.h>
  36 #include <TH1.h>
  37 #include <TVector2.h>
  38
  39 // Flow things
  40 #include "AliFlowEvent.h"
  41 #include "AliFlowTrack.h"
  42 #include "AliFlowV0.h"
  43 #include "AliFlowConstants.h"
  44 #include "AliFlowSelection.h"
  45 #include "AliFlowWeighter.h"
  46
  47 // ANSI things
  48 #include <math.h>
  49 #include <stdio.h>
  50 #include <stdlib.h>
  51 #include <iostream>
  52 using namespace std; //required for resolving the 'cout' symbol
  53
  54 ClassImp(AliFlowWeighter)
  55 //-----------------------------------------------------------------------
  56 AliFlowWeighter::AliFlowWeighter(const AliFlowSelection* flowSelect):
  57   fEventNumber(0), fTrackNumber(0), fNumberOfV0s(0), fNumberOfTracks(0), fPhiBins(0), fPhiMin(0.), fPhiMax(0.),
  58   fFlowEvent(0x0), fFlowTrack(0x0), fFlowSelect(0x0), fFlowTracks(0x0),
  59   fWgtFile(0x0), fWgtFileName("flowPhiWgt.hist.root"), fPhiWgtHistList(0x0)
  60 {
  61  // default constructor (selection given or default selection)
  62
  63  if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
  64  else           { fFlowSelect = new AliFlowSelection() ; }
  65
  66  // output file (histograms)
  67  //fWgtFile     = 0 ;
  68 }
  69 //-----------------------------------------------------------------------
  70 AliFlowWeighter::~AliFlowWeighter()
  71 {
  72  // default distructor (no actions)
  73 }
  74 //-----------------------------------------------------------------------
  75 Bool_t AliFlowWeighter::Init()
  76 {
  77 // sets some defaults for the analysis
  78
  79  cout << "* FlowWeighter *  -  Init()" << endl ; cout << endl ;
  80
  81  // Open output files (->plots)
  82  fWgtFile = new TFile(fWgtFileName.Data(), "RECREATE");
  83  fWgtFile->cd() ;
  84
  85  // counters and pointers to 0
  86  fEventNumber = 0 ;
  87  fTrackNumber = 0 ;
  88  fNumberOfTracks = 0 ;
  89  fNumberOfV0s    = 0 ;
  90  fFlowEvent  = 0 ;
  91  fFlowTrack  = 0 ;
  92  fFlowTracks = 0 ;
  93  //for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
  94
  95  // Histogram settings
  96  fPhiBins = AliFlowConstants::kPhiBins ;
  97  fPhiMin  = 0.;
  98  fPhiMax  = 2*TMath::Pi() ;
  99
 100  TString* histTitle ;
 101  for(int k = 0; k < AliFlowConstants::kSels; k++)
 102  {
 103   histTitle = new TString("Flow_BayPidMult_Sel");
 104   *histTitle += k+1;
 105   fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),AliFlowConstants::kPid,-0.5,((Float_t)AliFlowConstants::kPid-0.5));
 106   fHistFull[k].fHistBayPidMult->Sumw2() ;
 107   fHistFull[k].fHistBayPidMult->SetXTitle("e+/-  ,  mu+/-  ,  pi+/-  ,  K+/-  ,  p+/-  ,  d+/- ");
 108   fHistFull[k].fHistBayPidMult->SetYTitle("Counts");
 109   delete histTitle;
 110
 111   for(int j = 0; j < AliFlowConstants::kHars; j++)
 112   {
 113    // Tpc - Phi lab
 114    histTitle = new TString("Flow_Phi_TPC_Sel");
 115    *histTitle += k+1;
 116    histTitle->Append("_Har");
 117    *histTitle += j+1;
 118    fHistFull[k].fHistFullHar[j].fHistPhi = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 119    fHistFull[k].fHistFullHar[j].fHistPhi->SetXTitle("Azimuthal Angles (rad)");
 120    fHistFull[k].fHistFullHar[j].fHistPhi->SetYTitle("Counts");
 121    delete histTitle;
 122    // Tpc (plus)
 123    histTitle = new TString("Flow_Phi_TPCplus_Sel");
 124    *histTitle += k+1;
 125    histTitle->Append("_Har");
 126    *histTitle += j+1;
 127    fHistFull[k].fHistFullHar[j].fHistPhiPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 128    fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetXTitle("Azimuthal Angles (rad)");
 129    fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetYTitle("Counts");
 130    delete histTitle;
 131    // Tpc (minus)
 132    histTitle = new TString("Flow_Phi_TPCminus_Sel");
 133    *histTitle += k+1;
 134    histTitle->Append("_Har");
 135    *histTitle += j+1;
 136    fHistFull[k].fHistFullHar[j].fHistPhiMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 137    fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetXTitle("Azimuthal Angles (rad)");
 138    fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetYTitle("Counts");
 139    delete histTitle;
 140    // Tpc (cross)
 141    histTitle = new TString("Flow_Phi_TPCcross_Sel");
 142    *histTitle += k+1;
 143    histTitle->Append("_Har");
 144    *histTitle += j+1;
 145    fHistFull[k].fHistFullHar[j].fHistPhiAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 146    fHistFull[k].fHistFullHar[j].fHistPhiAll->SetXTitle("Azimuthal Angles (rad)");
 147    fHistFull[k].fHistFullHar[j].fHistPhiAll->SetYTitle("Counts");
 148    delete histTitle;
 149
 150    // Tpc - Phi lab flattened
 151    histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
 152    *histTitle += k+1;
 153    histTitle->Append("_Har");
 154    *histTitle += j+1;
 155    fHistFull[k].fHistFullHar[j].fHistPhiFlat = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 156    fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetXTitle("Azimuthal Angles (rad)");
 157    fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetYTitle("Counts");
 158    delete histTitle;
 159    // Tpc Plus
 160    histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
 161    *histTitle += k+1;
 162    histTitle->Append("_Har");
 163    *histTitle += j+1;
 164    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 165    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetXTitle("Azimuthal Angles (rad)");
 166    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetYTitle("Counts");
 167    delete histTitle;
 168    // Tpc Minus
 169    histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
 170    *histTitle += k+1;
 171    histTitle->Append("_Har");
 172    *histTitle += j+1;
 173    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 174    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetXTitle("Azimuthal Angles (rad)");
 175    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetYTitle("Counts");
 176    delete histTitle;
 177    // Tpc cross
 178    histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
 179    *histTitle += k+1;
 180    histTitle->Append("_Har");
 181    *histTitle += j+1;
 182    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 183    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetXTitle("Azimuthal Angles (rad)");
 184    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetYTitle("Counts");
 185    delete histTitle;
 186   }
 187  }
 188
 189  return kTRUE ;
 190 }
 191 //-----------------------------------------------------------------------
 192
 193 //-----------------------------------------------------------------------
 194 Bool_t AliFlowWeighter::Finish()
 195 {
 196  // Calls the method to fill wgt histograms, then saves them
 197  // on the fWgtFile and closes the file .
 198
 199  cout << "* FlowWeighter *  -  Finish()" << endl ; cout << endl ;
 200
 201  Weightening() ;
 202  fWgtFile->cd() ;
 203
 204  // Write PhiWgt histograms
 205  fPhiWgtHistList->Write();
 206  delete fPhiWgtHistList ;
 207
 208  // Write Bayesian Weights for P.Id.
 209  for(int k=0;k<AliFlowConstants::kSels;k++) { fHistFull[k].fHistBayPidMult->Write() ; }
 210
 211  // Write the AliFlowSelection object
 212  fFlowSelect->Write();
 213  delete fFlowSelect ;
 214
 215  fWgtFile->Close();
 216
 217  cout << "    Finish()  -  Wgt file closed : " << fWgtFileName.Data() << endl ; cout << endl ;
 218
 219  return kTRUE ;
 220 }
 221 //-----------------------------------------------------------------------
 222 // ###
 223 //-----------------------------------------------------------------------
 224 Bool_t AliFlowWeighter::WeightEvent(AliFlowEvent* fFlowEvent)
 225 {
 226  // Reads the AliFlowEvent (* fFlowEvent) and loops over the
 227  // AliFlowTraks to fill phi histograms .
 228
 229  cout << " AliFlowWeighter::WeightEvent(" << fFlowEvent << " )   -   " << fEventNumber << endl ;
 230  if(!fFlowEvent) { return kFALSE ; }
 231
 232  if(fFlowSelect->Select(fFlowEvent))     // event selected - here below the ANALYSIS FLAGS are setted -
 233  {
 234   fFlowTracks = fFlowEvent->TrackCollection() ;
 235   fNumberOfTracks = fFlowTracks->GetEntries() ;
 236   cout << "       event ID = " << fFlowEvent->EventID() << " :  found " << fNumberOfTracks << " AliFlowTracks . " << endl ;
 237   fFlowEvent->SetNoWgt() ;
 238   fFlowEvent->SetSelections(fFlowSelect) ; // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
 239
 240   TracksLoop(fFlowTracks) ;
 241  }
 242  else
 243  {
 244   cout << "       * " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ;
 245   delete fFlowEvent  ; fFlowEvent = 0 ;
 246   return kFALSE ;
 247  }
 248  fEventNumber++ ;
 249
 250  return kTRUE ;
 251 }
 252 //-----------------------------------------------------------------------
 253 // ###
 254 //-----------------------------------------------------------------------
 255 void AliFlowWeighter::TracksLoop(TObjArray* fFlowTracks)
 256 {
 257  // fills phi and PId histograms
 258
 259  cout << " Tracks Loop . " << endl ;
 260
 261  Float_t phi, eta, zFirstPoint ; // , zLastPoint ;
 262  Int_t  fPidId ;
 263  Char_t pid[10] = "0" ;
 264  for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++)
 265  {
 266   fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
 267   // cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ;
 268
 269   phi = fFlowTrack->Phi() ;
 270   eta = fFlowTrack->Eta() ;
 271   zFirstPoint = fFlowTrack->ZFirstPoint() ;
 272   //zLastPoint = fFlowTrack->ZLastPoint() ;
 273   strcpy(pid,fFlowTrack->Pid()) ;
 274
 275   fPidId = -1 ;
 276   if(strstr(pid,"e"))       { fPidId = 0 ; }
 277   else if(strstr(pid,"mu")) { fPidId = 1 ; }
 278   else if(strstr(pid,"pi")) { fPidId = 2 ; }
 279   else if(strstr(pid,"k"))  { fPidId = 3 ; }
 280   else if(strstr(pid,"pr")) { fPidId = 4 ; }
 281   else if(strstr(pid,"d"))  { fPidId = 5 ; }
 282
 283   // Looping over Selections and Harmonics
 284   for (int k = 0; k < AliFlowConstants::kSels; k++)
 285   {
 286    fFlowSelect->SetSelection(k) ;
 287    for (int j = 0; j < AliFlowConstants::kHars; j++)
 288    {
 289     fFlowSelect->SetHarmonic(j);
 290     if(fFlowSelect->Select(fFlowTrack))
 291     {
 292      Bool_t kTpcPlus  = kFALSE ;
 293      Bool_t kTpcMinus = kFALSE ;
 294      Bool_t kTpcAll   = kFALSE ;
 295
 296      // Set Tpc (+ and -)
 297      if(fFlowTrack->FitPtsTPC())
 298      {
 299       if(zFirstPoint >= 0. && eta > 0.)      { kTpcPlus  = kTRUE ; }
 300       else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
 301       else                                   { kTpcAll   = kTRUE ; }
 302      }
 303
 304      // Phi distribution (particle for R.P.)
 305      Float_t wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ;
 306      if(kTpcPlus)       { fHistFull[k].fHistFullHar[j].fHistPhiPlus->Fill(phi,wt) ;  }
 307      else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiMinus->Fill(phi,wt) ; }
 308      else if(kTpcAll)   { fHistFull[k].fHistFullHar[j].fHistPhiAll->Fill(phi,wt) ;   }
 309      fHistFull[k].fHistFullHar[j].fHistPhi->Fill(phi,wt) ;
 310
 311      // PID Multiplicities (particle for R.P.) - just once for each selection
 312      if(j==0)           { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
 313     }
 314    }
 315   }
 316  }
 317
 318  return ;
 319 }
 320 //-----------------------------------------------------------------------
 321 Bool_t AliFlowWeighter::Weightening()
 322 {
 323  // Calculates weights, and fills PhiWgt histograms .
 324  // This is called at the end of the event loop .
 325
 326  cout << " AliFlowWeighter::Weightening() " << endl ; cout << endl ;
 327
 328  // PhiWgt histogram collection
 329  fPhiWgtHistList = new TOrdCollection(4*AliFlowConstants::kSels*AliFlowConstants::kHars) ;
 330
 331  // Creates PhiWgt Histograms
 332  TString* histTitle ;
 333  for(Int_t k = 0; k < AliFlowConstants::kSels; k++)
 334  {
 335   for(Int_t j = 0; j < AliFlowConstants::kHars; j++)
 336   {
 337    // Tpc plus
 338    histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
 339    *histTitle += k+1;
 340    histTitle->Append("_Har");
 341    *histTitle += j+1;
 342    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 343    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Sumw2();
 344    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetXTitle("Azimuthal Angles (rad)");
 345    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetYTitle("PhiWgt");
 346    delete histTitle;
 347    // Tpc minus
 348    histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
 349    *histTitle += k+1;
 350    histTitle->Append("_Har");
 351    *histTitle += j+1;
 352    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 353    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Sumw2();
 354    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetXTitle("Azimuthal Angles (rad)");
 355    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetYTitle("PhiWgt");
 356    delete histTitle;
 357    // Tpc cross
 358    histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
 359    *histTitle += k+1;
 360    histTitle->Append("_Har");
 361    *histTitle += j+1;
 362    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 363    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Sumw2();
 364    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetXTitle("Azimuthal Angles (rad)");
 365    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetYTitle("PhiWgt");
 366    delete histTitle;
 367    // Tpc
 368    histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
 369    *histTitle += k+1;
 370    histTitle->Append("_Har");
 371    *histTitle += j+1;
 372    fHistFull[k].fHistFullHar[j].fHistPhiWgt = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
 373    fHistFull[k].fHistFullHar[j].fHistPhiWgt->Sumw2();
 374    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetXTitle("Azimuthal Angles (rad)");
 375    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetYTitle("PhiWgt");
 376    delete histTitle;
 377
 378    // Calculate PhiWgt
 379    Double_t meanPlus  = fHistFull[k].fHistFullHar[j].fHistPhiPlus->Integral() / (Double_t)fPhiBins ;
 380    Double_t meanMinus = fHistFull[k].fHistFullHar[j].fHistPhiMinus->Integral() / (Double_t)fPhiBins ;
 381    Double_t meanCross = fHistFull[k].fHistFullHar[j].fHistPhiAll->Integral() / (Double_t)fPhiBins ;
 382    Double_t meanTPC = fHistFull[k].fHistFullHar[j].fHistPhi->Integral() / (Double_t)fPhiBins ;
 383
 384    // Tpc
 385    for(Int_t i=0;i<fPhiBins;i++)
 386    {
 387     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinContent(i+1,meanPlus);
 388     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinError(i+1, 0.);
 389     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinContent(i+1,meanMinus);
 390     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinError(i+1, 0.);
 391     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinContent(i+1,meanCross);
 392     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinError(i+1, 0.);
 393     fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinContent(i+1,meanTPC);
 394     fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinError(i+1, 0.);
 395    }
 396
 397    if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " :  empty histogram ! " << endl ; }
 398    else
 399    {
 400     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiPlus);
 401     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiMinus);
 402     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Divide(fHistFull[k].fHistFullHar[j].fHistPhiAll);
 403     fHistFull[k].fHistFullHar[j].fHistPhiWgt->Divide(fHistFull[k].fHistFullHar[j].fHistPhi);
 404    }
 405
 406    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus);
 407    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus);
 408    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtAll);
 409    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgt);
 410   }
 411  }
 412
 413  return kTRUE ;
 414 }
 415 //-----------------------------------------------------------------------
 416 void AliFlowWeighter::PrintBayesian(Int_t selN) const
 417 {
 418  // Prints the normalized particle abundance of all events (selection selN).
 419  // Call this at the end of the loop, just before Finish() .
 420
 421  if(selN>AliFlowConstants::kSels) { selN = 0 ; }
 422  Char_t* names[AliFlowConstants::kPid] = {"e","mu","pi","k","p","d"} ;
 423  Double_t bayes = 0. ;
 424  Double_t totCount = (fHistFull[selN].fHistBayPidMult)->GetSumOfWeights() ;
 425  if(totCount)
 426  {
 427   cout << " Sel." << selN << " particles normalized abundance (tot. " << totCount << " tracks) : " ;
 428   for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++)
 429   {
 430    bayes = (fHistFull[selN].fHistBayPidMult->GetBinContent(ii+1) / totCount) ;
 431    cout << bayes << "_" << names[ii] << " ; " ;
 432   }
 433  }
 434  else {  cout << " Sel." << selN << " :  empty P.Id. histogram ! " << endl ; }
 435  cout << endl ;
 436
 437  return ;
 438 }
 439 //-----------------------------------------------------------------------
 440
 441
 442 #endif