PWG2/FLOW/AliCumulantsFunctions.cxx

   1 //*******************************
   2 // flow analysis with cumulants *
   3 // author: Ante Bilandzic       *
   4 // email: anteb@nikhef.nl       *
   5 //*******************************
   6
   7 #define AliCumulantsFunctions_cxx
   8
   9 #include "Riostream.h"
  10 #include "AliFlowCommonConstants.h"
  11 #include "TChain.h"
  12 #include "TFile.h"
  13 #include "TList.h" //NEW
  14 #include "TParticle.h"
  15
  16 #include "TProfile.h"
  17 #include "TProfile2D.h"
  18 #include "TProfile3D.h"
  19
  20 #include "AliFlowEventSimple.h"
  21 #include "AliFlowTrackSimple.h"
  22 #include "AliFlowAnalysisWithCumulants.h"
  23 #include "AliFlowCumuConstants.h"
  24 #include "AliFlowCommonConstants.h"
  25 #include "AliCumulantsFunctions.h"
  26
  27
  28 ClassImp(AliCumulantsFunctions)
  29
  30 //________________________________________________________________________
  31
  32 AliCumulantsFunctions::AliCumulantsFunctions():
  33   fIG(NULL),
  34   fDGRe(NULL),
  35   fDGIm(NULL),
  36   fifr(NULL),
  37   fdfr2(NULL),
  38   fdfr4(NULL),
  39   fdfr6(NULL),
  40   fdfr8(NULL),
  41   fAvMult(0)
  42  {
  43    //default constructor
  44  }
  45
  46 AliCumulantsFunctions::~AliCumulantsFunctions(){
  47   //desctructor
  48 }
  49
  50 AliCumulantsFunctions::AliCumulantsFunctions(TProfile2D *IntGen, TProfile3D *DiffGenRe, TProfile3D *DiffGenIm, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, Double_t CvM):
  51   fIG(IntGen),
  52   fDGRe(DiffGenRe),
  53   fDGIm(DiffGenIm),
  54   fifr(ifr),
  55   fdfr2(dfr2),
  56   fdfr4(dfr4),
  57   fdfr6(dfr6),
  58   fdfr8(dfr8),
  59   fAvMult(CvM)
  60  {
  61    //custom constructor
  62  }
  63
  64 //___________________________________________________________________________
  65 void AliCumulantsFunctions::Calculate(){
  66  //calculate final flow estimates
  67
  68      static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;//needed for numerics
  69      static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;//needed for numerics
  70      static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;//integrated flow coefficient to be calculated
  71      static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl;//the multiple in p=m*n (diff. flow)
  72      static const Int_t fgknBins=100;//number of pt bins //to be improved
  73      Double_t fR0=AliFlowCumuConstants::fgR0;//needed for numerics
  74      Double_t fPtMax=AliFlowCommonConstants::GetPtMax();//maximum pt
  75      Double_t fPtMin=AliFlowCommonConstants::GetPtMin();//minimum pt
  76      Double_t fBinWidth=(fPtMax-fPtMin)/fgknBins;//width of pt bin (in GeV)
  77
  78   Double_t fBvG[fgkPmax][fgkQmax]={0.};
  79
  80   for(Int_t p=0;p<fgkPmax;p++){
  81    for(Int_t q=0;q<fgkQmax;q++){
  82     fBvG[p][q]=fIG->GetBinContent(p+1,q+1);
  83    }
  84   }
  85
  86
  87   /////////////////////////////////////////////////////////////////////////////
  88   //////////////////gen. function for the cumulants////////////////////////////
  89   /////////////////////////////////////////////////////////////////////////////
  90
  91   Double_t fC[fgkPmax][fgkQmax]={0.};
  92   for (Int_t p=0;p<fgkPmax;p++){
  93     for (Int_t q=0;q<fgkQmax;q++){
  94       fC[p][q]=1.*fAvMult*(pow(fBvG[p][q],(1./fAvMult))-1.); //to be improved
  95     }
  96   }
  97
  98   /////////////////////////////////////////////////////////////////////////////
  99   ///////avaraging the gen. function for the cumulants over azimuth////////////
 100   /////////////////////////////////////////////////////////////////////////////
 101
 102   Double_t fCAv[fgkPmax]={0.};
 103   for (Int_t p=0;p<fgkPmax;p++){
 104     Double_t fTempHere=0.;
 105     for (Int_t q=0;q<fgkQmax;q++){
 106       fTempHere+=1.*fC[p][q];
 107     }
 108     fCAv[p]=1.*fTempHere/fgkQmax;
 109   }
 110
 111   /////////////////////////////////////////////////////////////////////////////
 112   //////////////////////////////////final results//////////////////////////////
 113   /////////////////////////////////////////////////////////////////////////////
 114
 115   Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
 116
 117   fCumulant[0]=(1./(fR0*fR0)) * (8.*fCAv[0] - 14.*fCAv[1] + (56./3.)*fCAv[2] - (35./2.)*fCAv[3] +
 118                               (56./5.)*fCAv[4] - (14./3.)*fCAv[5] + (8./7.)*fCAv[6] - (1./8.)*fCAv[7]);
 119
 120   fCumulant[1]=(1./pow(fR0,4.)) * ((-1924./35.)*fCAv[0] + (621./5.)*fCAv[1] - (8012./45.)*fCAv[2] +
 121                                  (691./4.)*fCAv[3] - (564./5.)*fCAv[4] + (2143./45.)*fCAv[5] -
 122                                  (412./35.)*fCAv[6] + (363./280.)*fCAv[7]);
 123
 124   fCumulant[2]=(1./pow(fR0,6.)) * (349.*fCAv[0] - (18353./20.)*fCAv[1] + (7173./5.)*fCAv[2] -
 125                                  1457.*fCAv[3] + (4891./5.)*fCAv[4] - (1683./4.)*fCAv[5] +
 126                                  (527./5.)*fCAv[6] - (469./40.)*fCAv[7]);
 127
 128   fCumulant[3]=(1./pow(fR0,8.)) * ((-10528./5.)*fCAv[0] + (30578./5.)*fCAv[1] - (51456./5.)*fCAv[2] +
 129                                  10993.*fCAv[3] - (38176./5.)*fCAv[4] + (16818./5.)*fCAv[5] -
 130                                  (4288./5.)*fCAv[6] + (967./10.)*fCAv[7]);
 131
 132   fCumulant[4]=(1./pow(fR0,10.)) * (11500.*fCAv[0] - 35800.*fCAv[1] + 63900.*fCAv[2] - 71600.*fCAv[3] +
 133                                   51620.*fCAv[4] - 23400.*fCAv[5] + 6100.*fCAv[6] - 700.*fCAv[7]);
 134
 135   fCumulant[5]=(1./pow(fR0,12.)) * (-52560.*fCAv[0] + 172080.*fCAv[1] - 321840.*fCAv[2] + 376200.*fCAv[3] -
 136                                   281520.*fCAv[4] + 131760.*fCAv[5] - 35280.*fCAv[6] + 4140.*fCAv[7]);
 137
 138   fCumulant[6]=(1./pow(fR0,14.)) * (176400.*fCAv[0] - 599760.*fCAv[1] + 1164240.*fCAv[2] - 1411200.*fCAv[3] +
 139                                   1093680.*fCAv[4] - 529200.*fCAv[5] + 146160.*fCAv[6] - 17640.*fCAv[7]);
 140
 141   fCumulant[7]=(1./pow(fR0,16.)) * (-322560*fCAv[0] + 1128960.*fCAv[1] - 2257920.*fCAv[2] + 2822400.*fCAv[3] -
 142                                   2257920.*fCAv[4] + 1128960.*fCAv[5] - 322560.*fCAv[6] + 40320.*fCAv[7]);
 143
 144
 145   cout<<""<<endl;
 146   cout<<"***************************"<<endl;
 147   cout<<"cumulants:"<<endl;
 148
 149   cout<<" c_"<<fgkFlow<<"{2} = "<<fCumulant[0]<<endl;
 150   cout<<" c_"<<fgkFlow<<"{4} = "<<fCumulant[1]<<endl;
 151   cout<<" c_"<<fgkFlow<<"{6} = "<<fCumulant[2]<<endl;
 152   cout<<" c_"<<fgkFlow<<"{8} = "<<fCumulant[3]<<endl;
 153   cout<<"c_"<<fgkFlow<<"{10} = "<<fCumulant[4]<<endl;
 154   cout<<"c_"<<fgkFlow<<"{12} = "<<fCumulant[5]<<endl;
 155   cout<<"c_"<<fgkFlow<<"{14} = "<<fCumulant[6]<<endl;
 156   cout<<"c_"<<fgkFlow<<"{16} = "<<fCumulant[7]<<endl;
 157
 158   cout<<""<<endl;
 159   cout<<"integrated flow: "<<endl;
 160
 161
 162   Double_t fV2=0.,fV4=0.,fV6=0.,fV8=0.;
 163   Double_t fSdQ[4]={0.};
 164   Double_t fChiQ[4]={0.};
 165    if (fCumulant[0]>=0.){
 166     fV2=sqrt(fCumulant[0]);
 167     //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
 168      //fChiQ[0]=fAvM*fV2/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.),0.5);
 169      //fSdQ[0]=pow(((1./(2.*fAvM*nEvents))*((1.+1.*pow(fChiQ[0],2))/(1.*pow(fChiQ[0],2)))),0.5);
 170      cout<<" v_"<<fgkFlow<<"{2} = "<<100.*fV2<<"%, chi{2} = "<<fChiQ[0]<<", sd{2} = "<<100.*fSdQ[0]<<"%"<<endl;//to be improved (2->fgkFlow)
 171      //fCommonHistsRes2->FillIntegratedFlow(100.*fV2,100.*fSdQ[0]);
 172      //fCommonHistsRes2->FillChi(fChiQ[0]);
 173      fifr->SetBinContent(1,100.*fV2);
 174     //}
 175    //} else {
 176     //cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
 177    } else {
 178     //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
 179    }
 180   if (fCumulant[1]<=0.){
 181     fV4=pow(-fCumulant[1],(1./4.));
 182     //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
 183      //fChiQ[1]=fAvM*fV4/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.),0.5);
 184      //fSdQ[1]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((1.+2.*pow(fChiQ[1],2)+(1./4.)*pow(fChiQ[1],4.)+(1./4.)*pow(fChiQ[1],6.))/((1./4.)*pow(fChiQ[1],6.)),.5);
 185      cout<<" v_"<<fgkFlow<<"{4} = "<<100.*fV4<<"%, chi{4} = "<<fChiQ[1]<<", sd{4} = "<<100.*fSdQ[1]<<"%"<<endl;//to be improved (2->fgkFlow)
 186      //fCommonHistsRes4->FillInteg8ratedFlow(100.*fV4,100.*fSdQ[1]);
 187      //fCommonHistsRes4->FillChi(fChiQ[1]);
 188      fifr->SetBinContent(2,100.*fV4);
 189     //} else {
 190       //cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
 191     //}
 192  // } else {
 193    // cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
 194   }
 195   if (fCumulant[2]>=0.){
 196     fV6=pow((1./4.)*fCumulant[2],(1./6.));
 197     //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.)>0.){
 198      //fChiQ[2]=fAvM*fV6/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.),0.5);
 199      //fSdQ[2]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((3.+18.*pow(fChiQ[2],2)+9.*pow(fChiQ[2],4.)+28.*pow(fChiQ[2],6.)+12.*pow(fChiQ[2],8.)+24.*pow(fChiQ[2],10.))/(24.*pow(fChiQ[2],10.)),.5);
 200      cout<<" v_"<<fgkFlow<<"{6} = "<<100.*fV6<<"%, chi{6} = "<<fChiQ[2]<<", sd{6} = "<<100.*fSdQ[2]<<"%"<<endl;//to be improved (2->fgkFlow)
 201      //fCommonHistsRes6->FillIntegratedFlow(100.*fV6,100.*fSdQ[2]);
 202      //fCommonHistsRes6->FillChi(fChiQ[2]);
 203      fifr->SetBinContent(3,100.*fV6);
 204     //} else {
 205      // cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
 206    // }
 207   } else {
 208     //cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
 209   }
 210   if (fCumulant[3]<=0.){
 211     fV8=pow(-(1./33.)*fCumulant[3],(1./8.));
 212     //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.)>0.){
 213      //fChiQ[3]=fAvM*fV8/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.),0.5);
 214      //fSdQ[3]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((12.+96.*pow(fChiQ[3],2)+72.*pow(fChiQ[3],4.)+304.*pow(fChiQ[3],6.)+257.*pow(fChiQ[3],8.)+804.*pow(fChiQ[3],10.)+363.*pow(fChiQ[3],12.)+726.*pow(fChiQ[3],14.))/(726.*pow(fChiQ[3],14.)),.5);
 215       cout<<" v_"<<fgkFlow<<"{8} = "<<100.*fV8<<"%, chi{8} = "<<fChiQ[3]<<", sd{8} = "<<100.*fSdQ[3]<<"%"<<endl;//to be improved (2->fgkFlow)
 216       //fCommonHistsRes8->FillIntegratedFlow(100.*fV8,100.*fSdQ[3]);
 217       //fCommonHistsRes8->FillChi(fChiQ[3]);
 218       fifr->SetBinContent(4,100.*fV8);
 219      //} else {
 220        //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;a
 221      //}
 222   } else {
 223     //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
 224   }
 225   if (fCumulant[4]>=0.){//fHistProIntFlow
 226     cout<<"v_"<<fgkFlow<<"{10} = "<<100.*pow((1./456.)*fCumulant[4],(1./10.))<<"%"<<endl;//to be improved (2->fgkFlow)
 227     fifr->SetBinContent(5,100.*pow((1./456.)*fCumulant[4],(1./10.)));
 228   } else {
 229       cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;  //to be improved (2->fgkFlow)
 230   }
 231   if (fCumulant[5]<=0.){
 232     cout<<"v_"<<fgkFlow<<"{12} = "<<100.*pow(-(1./9460.)*fCumulant[5],(1./12.))<<"%"<<endl;//to be improved (2->fgkFlow)
 233     fifr->SetBinContent(6,100.*pow(-(1./9460.)*fCumulant[5],(1./12.)));
 234   } else {
 235     cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl;  //to be improved (2->fgkFlow)
 236   }
 237   if (fCumulant[6]>=0.){
 238     cout<<"v_"<<fgkFlow<<"{14} = "<<100.*pow((1./274800.)*fCumulant[6],(1./14.))<<"%"<<endl;//to be improved (2->fgkFlow)
 239     fifr->SetBinContent(7,100.*pow((1./274800.)*fCumulant[6],(1./14.)));
 240   } else {
 241     cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;  //to be improved (2->fgkFlow)
 242   }
 243   if (fCumulant[7]<=0.){
 244     cout<<"v_"<<fgkFlow<<"{16} = "<<100.*pow(-(1./10643745.)*fCumulant[7],(1./16.))<<"%"<<endl;//to be improved (2->fgkFlow)
 245     fifr->SetBinContent(8,100.*pow(-(1./10643745.)*fCumulant[7],(1./16.)));
 246   } else {
 247     cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;  //to be improved (2->fgkFlow)
 248   }
 249   //cout<<"***************************"<<endl;
 250
 251   //DIFFERENTIAL FLOW CALCULATIONS STARTS HERE!!!
 252
 253   Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
 254   Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
 255
 256   for(Int_t b=0;b<fgknBins;b++){
 257     for(Int_t p=0;p<fgkPmax;p++){
 258       for(Int_t q=0;q<fgkQmax;q++){
 259         fX[b][p][q]=fDGRe->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
 260         fY[b][p][q]=fDGIm->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
 261       }
 262     }
 263   }
 264
 265
 266   Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
 267
 268   for (Int_t b=0;b<fgknBins;b++){
 269     for (Int_t p=0;p<fgkPmax;p++){
 270       Double_t fTempHere3=0.;
 271       for (Int_t q=0;q<fgkQmax;q++){
 272         fTempHere3+=cos(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fX[b][p][q] + sin(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fY[b][p][q];
 273       }
 274       fD[b][p]=1.*(pow(fR0*pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
 275       if(fD[b][p]){
 276        //cout<<"And this "<<b<<" "<<fD[b][p]<<endl;
 277       }
 278     }
 279   }
 280
 281   Double_t fDiffCumulant2[fgknBins]={0.};//implementing relation (12) from PG
 282   Double_t fDiffCumulant4[fgknBins]={0.};
 283   Double_t fDiffCumulant6[fgknBins]={0.};
 284   Double_t fDiffCumulant8[fgknBins]={0.};
 285
 286   for (Int_t b=0;b<fgknBins;b++){
 287     fDiffCumulant2[b]=(1./(fR0*fR0))*(4.*fD[b][0]-3.*fD[b][1]+(4./3.)*fD[b][2]-(1./4.)*fD[b][3]);
 288     fDiffCumulant4[b]=(1./pow(fR0,4.))*((-26./3.)*fD[b][0]+(19./2.)*fD[b][1]-(14./3.)*fD[b][2]+(11./12.)*fD[b][3]);
 289     fDiffCumulant6[b]=(1./pow(fR0,6.))*(18.*fD[b][0]-24.*fD[b][1]+14.*fD[b][2]-3.*fD[b][3]);
 290     fDiffCumulant8[b]=(1./pow(fR0,8.))*(-24.*fD[b][0]+36.*fD[b][1]-24.*fD[b][2]+6.*fD[b][3]);
 291   }
 292
 293   Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
 294   //Double_t fAvPt[fgknBins];
 295   Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];
 296
 297   cout<<"number of pt bins: "<<fgknBins<<endl;
 298   cout<<"****************************************"<<endl;
 299   for (Int_t b=0;b<fgknBins;b++){
 300     //if(fBinNoOfParticles[b]==0)continue;
 301     //fAvPt[b]=fBinMeanPt[b]/fBinNoOfParticles[b];
 302     cout<<"pt bin: "<<b*fBinWidth<<"-"<<(b+1)*fBinWidth<<" GeV"<<endl;
 303     //cout<<"number of particles in this pt bin: "<<tempNo->GetBinContent(b+1,1,1)<<endl;
 304     //cout<<"mean pt in this bin: "<<fAvPt[b]<<" GeV"<<endl;
 305     if(fCumulant[0]>=0){
 306       fv2[b]=100.*fDiffCumulant2[b]/pow(fCumulant[0],.5);
 307       //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
 308        //fSddiff2[b]=pow((1./(2.*fBinNoOfParticles[b]))*((1.+pow(fChiQ[0],2.))/pow(fChiQ[0],2.)),0.5);
 309        cout<<"v'_2/2{2} = "<<fv2[b]<<"%, "<<" "<<"sd{2} = "<<100.*fSddiff2[b]<<"%"<<endl;
 310        //fDiffFlowResults2->SetBinContent(b+1,fv2[b],100.*fSddiff2[b]);
 311        fdfr2->SetBinContent(b+1,fv2[b]);
 312       //} else {
 313          //cout<<"v'_2/2{2} = Im"<<endl;
 314       //}
 315     }else{
 316       cout<<"v'_2/2{2} = Im"<<endl;
 317     }
 318     if(fCumulant[1]<=0){
 319       fv4[b]=-100.*fDiffCumulant4[b]/pow(-fCumulant[1],.75);
 320       //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
 321        //fSddiff4[b]=pow((1./(2.*fBinNoOfParticles[b]))*((2.+6.*pow(fChiQ[1],2.)+pow(fChiQ[1],4.)+pow(fChiQ[1],6.))/pow(fChiQ[1],6.)),0.5);
 322        cout<<"v'_2/2{4} = "<<fv4[b]<<"%, "<<" "<<"sd{4} = "<<100.*fSddiff4[b]<<"%"<<endl;
 323        //fCommonHistsRes4->FillDifferentialFlow(b+1,fv4[b],100.*fSddiff4[b]);
 324        fdfr4->SetBinContent(b+1,fv4[b]);
 325       //} else {
 326        // cout<<"v'_2/2{4} = Im"<<endl;
 327       //}
 328     }else{
 329       cout<<"v'_2/2{4} = Im"<<endl;
 330     }
 331     if(fCumulant[2]>=0){
 332       cout<<"v'_2/2{6} = "<<100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)))<<"%"<<endl;
 333       fv6[b]=100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)));
 334       //fCommonHistsRes6->FillDifferentialFlow(b+1,fv6[b],0.);
 335       fdfr6->SetBinContent(b+1,fv6[b]);
 336     }else{
 337       cout<<"v'_2/2{6} = Im"<<endl;
 338     }
 339     if(fCumulant[3]<=0){
 340       cout<<"v'_2/2{8} = "<<-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.)))<<"%"<<endl;
 341       fv8[b]=-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.)));
 342       //fCommonHistsRes8->FillDifferentialFlow(b+1,fv8[b],0.);
 343       fdfr8->SetBinContent(b+1,fv8[b]);
 344     }else{
 345       cout<<"v'_2/2{8} = Im"<<endl;
 346     }
 347     cout<<"****************************************"<<endl;
 348   }
 349 }
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364