PWG2/FLOW/other/ShinIchi/AnaEveMyTree.C

   1 #define AnaEveMyTree_cxx
   2 #include "AnaEveMyTree.h"
   3 #include <TH2.h>
   4 #include <TStyle.h>
   5 #include <TCanvas.h>
   6
   7 void AnaEveMyTree::Loop()
   8 {
   9 //      Root > .L AnaEveMyTree.C
  10 //      Root > AnaEveMyTree t
  11 //      Root > t.Loop();       // Loop on all entries
  12 //
  13 // T0C -3.3 ~ -2.9
  14 // T0A  4.5 ~  5.0
  15 //
  16 // V0C -1.7 ~ -2.2 ~ -2.7 ~ -3.2 ~ -3.7
  17 // V0A  2.8 ~  3.4 ~  3.9 ~  4.5 ~  5.1
  18 //
  19  float pi=acos(-1.0);
  20  float t0phi[24];
  21  for(int it=0; it<12; it++) {
  22   t0phi[it]=(it+1.0)*2.0*pi/12.0;
  23   t0phi[it+12]=it*2.0*pi/12.0;
  24  }
  25  for(int it=0; it<24; it++) {
  26   float phi=t0phi[it];
  27   t0phi[it]=atan2(sin(phi),cos(phi));
  28  }
  29  float v0phi[64];
  30  int is=0;
  31  for(int iRing = 0; iRing<4 ; iRing++){
  32   for(int iSec = 0; iSec<8 ; iSec++){
  33    v0phi[is] = 22.5 + 45. * iSec;
  34    is++;
  35   }
  36  }
  37  for(int iRing = 0; iRing<4 ; iRing++){
  38   for(int iSec = 0; iSec<8 ; iSec++){
  39    v0phi[is] = 22.5 + 45. * iSec;
  40    is++;
  41   }
  42  }
  43  for (int i=0; i<64; i++) {
  44   float phi=v0phi[i]*pi/180.0;
  45   v0phi[i]=atan2(sin(phi),cos(phi));
  46   cout << v0phi[i] << " ";
  47   if (i%8==7) cout << endl;
  48  }
  49  float z0gai[8],t0gai[24],v0gai[64];
  50  float mean[10][10][4][27][2];
  51  float widt[10][10][4][27][2];
  52  float four[10][10][4][27][2][8];
  53  float f0,f1,f2,f3,f4,f5,f6,f7;
  54  ifstream ifs;
  55  ifs.open("AnaEveMyTree.cal");
  56  for (int i=0; i<1; i++) {
  57   ifs>>f0>>f1>>f2>>f3>>f4>>f5>>f6>>f7;
  58   z0gai[i*8+0]=f0; z0gai[i*8+1]=f1;
  59   z0gai[i*8+2]=f2; z0gai[i*8+3]=f3;
  60   z0gai[i*8+4]=f4; z0gai[i*8+5]=f5;
  61   z0gai[i*8+6]=f6; z0gai[i*8+7]=f7;
  62   cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
  63        << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
  64  }
  65  for (int i=0; i<3; i++) {
  66   ifs>>f0>>f1>>f2>>f3>>f4>>f5>>f6>>f7;
  67   t0gai[i*8+0]=f0; t0gai[i*8+1]=f1;
  68   t0gai[i*8+2]=f2; t0gai[i*8+3]=f3;
  69   t0gai[i*8+4]=f4; t0gai[i*8+5]=f5;
  70   t0gai[i*8+6]=f6; t0gai[i*8+7]=f7;
  71   cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
  72        << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
  73  }
  74  for (int i=0; i<8; i++) {
  75   ifs>>f0>>f1>>f2>>f3>>f4>>f5>>f6>>f7;
  76   v0gai[i*8+0]=f0; v0gai[i*8+1]=f1;
  77   v0gai[i*8+2]=f2; v0gai[i*8+3]=f3;
  78   v0gai[i*8+4]=f4; v0gai[i*8+5]=f5;
  79   v0gai[i*8+6]=f6; v0gai[i*8+7]=f7;
  80   cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
  81        << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
  82  }
  83  for (int ic=0; ic<10; ic++) {
  84   for (int iz=0; iz<10; iz++) {
  85    for (int ih=0; ih<4; ih++) {
  86     for (int id=0; id<27; id++) {
  87      for (int ib=0; ib<2; ib++) {
  88       mean[ic][iz][ih][id][ib]=0.0;
  89       widt[ic][iz][ih][id][ib]=1.0;
  90      }
  91      ifs >> f0 >> f1 >> f2 >> f3;
  92      if (f1==0.0 || f1<0.0) f1=1.0;
  93      if (f3==0.0 || f3<0.0) f3=1.0;
  94      mean[ic][iz][ih][id][0]=f0;
  95      widt[ic][iz][ih][id][0]=f1;
  96      mean[ic][iz][ih][id][1]=f2;
  97      widt[ic][iz][ih][id][1]=f3;
  98      if (ic==4 && iz==4 && ih==0) {
  99       cout << f0 << " " << f1 << " " << f2 << " " << f3 << endl;
 100      }
 101      for (int ib=0; ib<2; ib++) {
 102       for (int io=0; io<8; io++) {
 103        four[ic][iz][ih][id][ib][io]=0.0;
 104       }
 105       ifs >> f0 >> f1 >> f2 >> f3 >> f4 >> f5 >> f6 >> f7;
 106       four[ic][iz][ih][id][ib][0]=f0;
 107       four[ic][iz][ih][id][ib][1]=f1;
 108       four[ic][iz][ih][id][ib][2]=f2;
 109       four[ic][iz][ih][id][ib][3]=f3;
 110       four[ic][iz][ih][id][ib][4]=f4;
 111       four[ic][iz][ih][id][ib][5]=f5;
 112       four[ic][iz][ih][id][ib][6]=f6;
 113       four[ic][iz][ih][id][ib][7]=f7;
 114       if (ic==4 && iz==4 && ih==0) {
 115        cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
 116             << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
 117       }
 118      }
 119     }
 120    }
 121   }
 122  }
 123  ifs.close();
 124  float av1=0;
 125  float no1=0;
 126  for (int i=0; i<8;  i++) {
 127   if (z0gai[i]>3) {
 128    av1+=z0gai[i];
 129    no1+=1.0;
 130   } else {
 131    z0gai[i]=-1.0;
 132   }
 133  }
 134  av1/=no1;
 135  float av2=0;
 136  float no2=0;
 137  for (int i=0; i<24;  i++) {
 138   if (t0gai[i]>3) {
 139    av2+=t0gai[i];
 140    no2+=1.0;
 141   } else {
 142    t0gai[i]=-1.0;
 143   }
 144  }
 145  av2/=no2;
 146  float av3=0;
 147  float no3=0;
 148  for (int i=0; i<64;  i++) {
 149   if (v0gai[i]>3) {
 150    av3+=v0gai[i];
 151    no3+=1.0;
 152   } else {
 153    v0gai[i]=-1.0;
 154   }
 155  }
 156  av3/=no3;
 157  cout << av1 << " " << av2 << " " << av3 << endl;
 158  TFile *tf = new TFile("AnaEveMyTree.root","recreate");
 159  TH2D     *cnzv = new TH2D("cnzv","",50,0.0,11000.0,50,-15.0,15.0);
 160  TProfile *det0g = new TProfile("det0g","", 16,0.0, 16.0   ,-600.0,6000.0);
 161  TH2D     *det0m = new TH2D    ("det0m","", 16,0.0, 16.0,50,-300.0,3000.0);
 162  TProfile *det1g = new TProfile("det1g","", 48,0.0, 48.0   , -20.0, 200.0);
 163  TH2D     *det1m = new TH2D    ("det1m","", 48,0.0, 48.0,50, -10.0, 100.0);
 164  TProfile *det2g = new TProfile("det2g","",128,0.0,128.0   ,-100.0,1000.0);
 165  TH2D     *det2m = new TH2D    ("det2m","",128,0.0,128.0,50, -50.0, 500.0);
 166  char name[80];
 167  TH2D     *zdcxy[2][5];
 168  for (int i=0; i<2; i++) {
 169   for (int j=0; j<5; j++) {
 170    sprintf(name,"zdcxy%d%d",i,j);
 171    zdcxy[i][j] = new TH2D(name,"",50,-2,2,50,-2,2);
 172   }
 173  }
 174  float rng0=11000.0;
 175  float rng1=av3*64.0*4.5;
 176  float rng2=av2*12.0*7.0;
 177  float rng3=av1* 8.0*2.2;
 178  TH2D *cor00 = new TH2D("corr00","",50,0.0,rng0,50,0.0,rng1);
 179  TH2D *cor01 = new TH2D("corr01","",50,0.0,rng0,50,0.0,rng2);
 180  TH2D *cor02 = new TH2D("corr02","",50,0.0,rng0,50,0.0,rng3);
 181  TH2D *cor03 = new TH2D("corr03","",50,0.0,rng1,50,0.0,rng2);
 182  TH2D *cor04 = new TH2D("corr04","",50,0.0,rng1,50,0.0,rng3);
 183  TH2D *cor05 = new TH2D("corr05","",50,0.0,rng2,50,0.0,rng3);
 184  TH2D *cor10 = new TH2D("corr10","",50,0.0,rng0,50,0.0,rng1);
 185  TH2D *cor11 = new TH2D("corr11","",50,0.0,rng0,50,0.0,rng2);
 186  TH2D *cor12 = new TH2D("corr12","",50,0.0,rng0,50,0.0,rng3);
 187  TH2D *cor13 = new TH2D("corr13","",50,0.0,rng1,50,0.0,rng2);
 188  TH2D *cor14 = new TH2D("corr14","",50,0.0,rng1,50,0.0,rng3);
 189  TH2D *cor15 = new TH2D("corr15","",50,0.0,rng2,50,0.0,rng3);
 190  TProfile *ave[10][10][4][27];
 191  TProfile *flt[10][10][4][27];
 192  for (int ic=0; ic<10; ic++) {
 193   for (int iz=0; iz<10; iz++) {
 194    for (int ih=0; ih<4; ih++) {
 195     for (int id=0; id<27; id++) {
 196      sprintf(name,"ave%d%d%d%d",ic,iz,ih,id);
 197      ave[ic][iz][ih][id] = new TProfile(name,name,4,-0.5,3.5,-1.0E+06,1.0E+06,"S");
 198      sprintf(name,"flt%d%d%d%d",ic,iz,ih,id);
 199      flt[ic][iz][ih][id] = new TProfile(name,name,32,-0.5,31.5,-1.0E+06,1.0E+06);
 200     }
 201    }
 202   }
 203  }
 204  TProfile *res[4][22][2];
 205  TH2D     *cor[4][22][5];
 206  for (int ih=0; ih<4; ih++) {
 207   for (int id=0; id<22; id++) {
 208    for (int ia=0; ia<2; ia++) {
 209     sprintf(name,"res%d%d%d",ih,id,ia);
 210     res[ih][id][ia] = new TProfile(name,name,10,-0.5,9.5,-1.1,1.1);
 211    }
 212    for (int ia=0; ia<5; ia++) {
 213     sprintf(name,"cor%d%d%d",ih,id,ia);
 214     cor[ih][id][ia] = new TH2D(name,name,50,-pi,pi,50,-pi,pi);
 215    }
 216   }
 217  }
 218  TH1D *t0d[24][5],*t0e[24][5];
 219  for (int i=0; i<24; i++) {
 220   for (int j=0; j<5; j++) {
 221    if (j<2) {
 222     sprintf(name,"t0d%d_%d",i,j);
 223     t0d[i][j] = new TH1D(name,"",24,-pi,pi);
 224     sprintf(name,"t0e%d_%d",i,j);
 225     t0e[i][j] = new TH1D(name,"",24,-pi,pi);
 226    } else {
 227     sprintf(name,"t0d%d_%d",i,j);
 228     t0d[i][j] = new TH1D(name,"",24,-pi/2.,pi/2.);
 229     sprintf(name,"t0e%d_%d",i,j);
 230     t0e[i][j] = new TH1D(name,"",24,-pi/2.,pi/2.);
 231    }
 232   }
 233  }
 234  TH2D *tvcr0[10][4],*tvcr1[10][4];
 235  for (int i=0; i<10; i++) {
 236   for (int j=0; j<4; j++) {
 237    sprintf(name,"tvcr0_%d_%d",i,j);
 238    tvcr0[i][j] = new TH2D(name,"",48,0,48,24,-pi,pi);
 239    sprintf(name,"tvcr1_%d_%d",i,j);
 240    tvcr1[i][j] = new TH2D(name,"",128,0,128,24,-pi,pi);
 241   }
 242  }
 243  int neve[10][10];
 244  float buft[10][10][3][24];
 245  float bufv[10][10][3][64];
 246  for (int i=0; i<10; i++) {
 247   for (int j=0; j<10; j++) {
 248    neve[i][i]=0;
 249   }
 250  }
 251  if (fChain == 0) return;
 252  Long64_t nentries = fChain->GetEntriesFast();
 253  cout << nentries << endl;
 254  Long64_t nbytes = 0, nb = 0;
 255  int ievent==0;
 256  for (Long64_t jentry=0; jentry<nentries;jentry++) {
 257   Long64_t ientry = LoadTree(jentry);
 258   if (ientry < 0) break;
 259   nb = fChain->GetEntry(jentry);   nbytes += nb;
 260   if (t1s >   1 && t2s >  1
 261    && v1s >  10 && v2s > 10
 262    && z1s > 100 && z2s >100
 263    && zvt > -10 && zvt < 10
 264    && v2s < 0.8*ntr+1000
 265    && v2s > 0.8*ntr-1000
 266    && t1s > 0.02*ntr) {
 267   cnzv->Fill(ntr,zvt);
 268   int icen=(int)(ntr/900.0);
 269   if (icen<0) icen=0;
 270   if (icen>9) icen=9;
 271   int izvt=(int)((zvt+10.0)/2.0);
 272   if (izvt<0) izvt=0;
 273   if (izvt>9) izvt=9;
 274
 275   float sumxy[4][27][4];
 276   for (int ih=0; ih<4; ih++) {
 277    for (int id=0; id<27; id++) {
 278     for (int iv=0; iv<4; iv++) {
 279      sumxy[ih][id][iv]=0;
 280     }
 281    }
 282   }
 283   float z0co[8];
 284   float z0sm=0;
 285   float z0x[8] = {-1,  1, -1,  1,  1, -1,  1, -1};
 286   float z0y[8] = {-1, -1,  1,  1, -1, -1,  1,  1};
 287   for (int i=0; i<8; i++) {
 288    float gain=z0gai[i];
 289    if (gain<0) gain=0;
 290    else gain=av1/gain;
 291    if (i<4) { // for C-side
 292     det0g->Fill(i+0.5,z1nt[i+1]);
 293     det0m->Fill(i+0.5,z1nt[i+1]);
 294     det0g->Fill(i+8.5,z1nt[i+1]*gain);
 295     det0m->Fill(i+8.5,z1nt[i+1]*gain);
 296     z0co[i]=z1nt[i+1]*gain;
 297     z0sm+=z1nt[i+1]*gain;
 298    } else { // for A-side
 299     det0g->Fill(i+0.5,z2nt[i-3]);
 300     det0m->Fill(i+0.5,z2nt[i-3]);
 301     det0g->Fill(i+8.5,z2nt[i-3]*gain);
 302     det0m->Fill(i+8.5,z2nt[i-3]*gain);
 303     z0co[i]=z2nt[i-3]*gain;
 304     z0sm+=z2nt[i-3]*gain;
 305    }
 306    sumxy[0][i/4][0]+=z0co[i]*z0x[i];
 307    sumxy[0][i/4][1]+=z0co[i]*z0y[i];
 308    sumxy[0][i/4][2]+=z0co[i];
 309    float sgn=-1.0*(1-(int)(i/4)); // sign flip for C-side
 310    sumxy[0][2][0]+=z0co[i]*z0x[i]*sgn;
 311    sumxy[0][2][1]+=z0co[i]*z0y[i]*sgn;
 312    sumxy[0][2][2]+=z0co[i];
 313   }
 314   for (int i=0; i<3; i++) {
 315    sumxy[0][i][0]/=sumxy[0][i][2];
 316    sumxy[0][i][1]/=sumxy[0][i][2];
 317   }
 318   zdcxy[0][icen/2]->Fill(sumxy[0][0][0],sumxy[0][1][0]);
 319   zdcxy[1][icen/2]->Fill(sumxy[0][0][1],sumxy[0][1][1]);
 320   float t0co[24];
 321   float t0sm=0;
 322   for (int i=0; i<24; i++) {
 323    float gain=t0gai[i];
 324    if (gain<0) gain=0;
 325    else gain=av2/gain;
 326    if (t0tim[i]>100 && t0amp[i]>0) {
 327     det1g->Fill(i+0.5,t0amp[i]);
 328     det1m->Fill(i+0.5,t0amp[i]);
 329     det1g->Fill(i+24.5,t0amp[i]*gain);
 330     det1m->Fill(i+24.5,t0amp[i]*gain);
 331     t0co[i]=t0amp[i]*gain;
 332     t0sm+=t0amp[i]*gain;
 333    } else {
 334     t0co[i]=0.0;
 335    }
 336   }
 337   float v0co[64];
 338   float v0sm=0;
 339   for (int i=0; i<64; i++) {
 340    float gain=v0gai[i];
 341    if (gain<0) gain=0;
 342    else gain=av3/gain;
 343    if (v0mul[i]>0) {
 344     det2g->Fill(i+0.5,v0mul[i]);
 345     det2m->Fill(i+0.5,v0mul[i]);
 346     det2g->Fill(i+64.5,v0mul[i]*gain);
 347     det2m->Fill(i+64.5,v0mul[i]*gain);
 348     v0co[i]=v0mul[i]*gain;
 349     v0sm+=v0mul[i]*gain;
 350    } else {
 351     v0co[i]=0.0;
 352    }
 353   }
 354   for (int i=0; i<14; i++) {
 355    int id=(int)(i/12);
 356    for (int ih=0; ih<4; ih++) {
 357     float sgn=1;
 358     if (ih%2==0) sgn=-1.0*(1-id); // sign flip for C-side
 359     sumxy[ih][3+id][0]+=t0co[i]*cos((ih+1.0)*t0phi[i]);
 360     sumxy[ih][3+id][1]+=t0co[i]*sin((ih+1.0)*t0phi[i]);
 361     sumxy[ih][3+id][2]+=t0co[i];
 362     sumxy[ih][5][0]+=t0co[i]*cos((ih+1.0)*t0phi[i]);
 363     sumxy[ih][5][1]+=t0co[i]*sin((ih+1.0)*t0phi[i]);
 364     sumxy[ih][5][2]+=t0co[i];
 365     sumxy[ih][6][0]+=t0co[i]*cos((ih+1.0)*t0phi[i])*sgn;
 366     sumxy[ih][6][1]+=t0co[i]*sin((ih+1.0)*t0phi[i])*sgn;
 367     sumxy[ih][6][2]+=t0co[i];
 368    }
 369   }
 370   for (int i=0; i<64; i++) {
 371    int id=(int)(i/32);
 372    int jd=(int)(i/8);
 373    int kd=jd%4;
 374    for (int ih=0; ih<4; ih++) {
 375     float sgn=1;
 376     if (ih%2==0) sgn=-1.0*(1-id); // sign flip for C-side
 377     sumxy[ih][7+id][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
 378     sumxy[ih][7+id][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
 379     sumxy[ih][7+id][2]+=v0co[i];
 380     sumxy[ih][9][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
 381     sumxy[ih][9][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
 382     sumxy[ih][9][2]+=v0co[i];
 383     sumxy[ih][10][0]+=v0co[i]*cos((ih+1.0)*v0phi[i])*sgn;
 384     sumxy[ih][10][1]+=v0co[i]*sin((ih+1.0)*v0phi[i])*sgn;
 385     sumxy[ih][10][2]+=v0co[i];
 386     sumxy[ih][11+4*kd+id][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
 387     sumxy[ih][11+4*kd+id][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
 388     sumxy[ih][11+4*kd+id][2]+=v0co[i];
 389     sumxy[ih][13+4*kd][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
 390     sumxy[ih][13+4*kd][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
 391     sumxy[ih][13+4*kd][2]+=v0co[i];
 392     sumxy[ih][14+4*kd][0]+=v0co[i]*cos((ih+1.0)*v0phi[i])*sgn;
 393     sumxy[ih][14+4*kd][1]+=v0co[i]*sin((ih+1.0)*v0phi[i])*sgn;
 394     sumxy[ih][14+4*kd][2]+=v0co[i];
 395    }
 396   }
 397   if (icen>0 && icen<7) {
 398    for (int i=0; i<24; i++) {
 399     float amp1=t0co[i]/av2;
 400     float phi1=t0phi[i];
 401     for (int j=0; j<64; j++) {
 402      int id=(int)(j/32);
 403      float amp2=v0co[j]/av3;
 404      float phi2=v0phi[j];
 405      float dphi=phi2-phi1;
 406      dphi=atan2(sin(dphi),cos(dphi));
 407      t0d[i][id]->Fill(phi2,amp1*amp2);
 408      t0e[i][id]->Fill(dphi,amp1*amp2);
 409     }
 410    }
 411   }
 412   int ieve=neve[icen][izvt];
 413   for (int i=0; i<24; i++) buft[icen][izvt][ieve][i]=t0co[i]/av2;
 414   for (int i=0; i<64; i++) bufv[icen][izvt][ieve][i]=v0co[i]/av3;
 415   neve[icen][izvt]++;
 416   if (neve[icen][izvt]==3) {
 417    neve[icen][izvt]=0;
 418    for (int ie=0; ie<3; ie++) {
 419     for (int je=0; je<3; je++) {
 420      int mix=1;
 421      if (ie==je) mix=0;
 422      for (int ip=0; ip<24; ip++) {
 423       int iq=(int)(ip/12);
 424       float phi1=t0phi[ip];
 425       float amp1=buft[icen][izvt][ie][ip];
 426       for (int jp=0; jp<64; jp++) {
 427        int jq=(int)(jp/32);
 428        float phi2=v0phi[jp];
 429        float amp2=bufv[icen][izvt][je][jp];
 430        float dphi=phi2-phi1;
 431        dphi=atan2(sin(dphi),cos(dphi));
 432        tvcr0[icen][iq*2+jq]->Fill(mix*24+ip+0.5,dphi,amp1*amp2);
 433        tvcr1[icen][iq*2+jq]->Fill(mix*64+jp+0.5,dphi,amp1*amp2);
 434       }
 435      }
 436     }
 437    }
 438   }
 439
 440 // r.p. flattening // -----------------------------------
 441   int calFlag=1;
 442   for (int ih=0; ih<4; ih++) {
 443    for (int id=0; id<27; id++) {
 444     for (int ib=0; ib<2; ib++) {
 445      if (calFlag>0) ave[icen][izvt][ih][id]->Fill(ib+0.0,sumxy[ih][id][ib]);
 446      float sxy=sumxy[ih][id][ib];
 447      float mxy=mean[icen][izvt][ih][id][ib];
 448      float wxy=widt[icen][izvt][ih][id][ib];
 449      sumxy[ih][id][ib]=(sxy-mxy)/wxy;
 450      if (calFlag>0) ave[icen][izvt][ih][id]->Fill(ib+2.0,sumxy[ih][id][ib]);
 451     }
 452     if (sumxy[ih][id][2]>0.0) {
 453      sumxy[ih][id][3]=atan2(sumxy[ih][id][1],sumxy[ih][id][0])/(ih+1.0);
 454      float psi=sumxy[ih][id][3]*(ih+1.0);
 455      float dp=0.0;
 456      for (int io=0; io<8; io++) {
 457       float cc=cos((io+1.0)*psi);
 458       float ss=sin((io+1.0)*psi);
 459       if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+0.0,cc);
 460       if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+8.0,ss);
 461       float aa=four[icen][izvt][ih][id][0][io]; // mean cos
 462       float bb=four[icen][izvt][ih][id][1][io]; // mean sin
 463       dp+=(aa*ss-bb*cc)*2.0/(io+1.0);
 464      }
 465      psi+=dp;
 466      for (int io=0; io<8; io++) {
 467       float cc=cos((io+1.0)*psi);
 468       float ss=sin((io+1.0)*psi);
 469       if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+16.0,cc);
 470       if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+24.0,ss);
 471      }
 472      sumxy[ih][id][3]=atan2(sin(psi),cos(psi))/(ih+1.0);
 473     } else {
 474      sumxy[ih][id][3]=-9999.9;
 475     }
 476    }
 477   }
 478
 479   if (icen>0 && icen<7) {
 480    for (int i=0; i<24; i++) {
 481     float amp1=t0co[i]/av2;
 482     float phi1=t0phi[i];
 483     for (int id=7; id<10; id++) {
 484      float phi2=sumxy[1][id][3];
 485      if (phi2>-9000) {
 486       float dphi=phi2-phi1;
 487       dphi=atan2(sin(2.0*dphi),cos(2.0*dphi))/2.0;
 488       t0d[i][id-5]->Fill(phi2,amp1);
 489       t0e[i][id-5]->Fill(dphi,amp1);
 490      }
 491     }
 492    }
 493   }
 494
 495 // r.p. correlation // -----------------------------------
 496   for (int ih=0; ih<4; ih++) {
 497    for (int id=0; id<22; id++) {
 498     int id1=0;
 499     int id2=0;
 500     int ih1=ih;
 501     int ih2=ih;
 502     if (id==0) {id1=0; id2=1; ih1=0; ih2=0;} // z0c-z0c (1st)
 503     if (id==1) {id1=3; id2=4;}               // t0c-t0a
 504     if (id==2) {id1=7; id2=8;}               // v0c-v0a
 505     if (id==3) {id1=2; id2=5; ih1=0;}        // zd1-t0ac  6 for t0ac flip
 506     if (id==4) {id1=2; id2=9; ih1=0;}        // zd1-v0ac 10 for v0ac flip
 507     if (id==5) {id1=3; id2=7;}               // t0c-v0c
 508     if (id==6) {id1=3; id2=8;}               // t0c-v0a
 509     if (id==7) {id1=4; id2=7;}               // t0a-v0c
 510     if (id==8) {id1=4; id2=8;}               // t0a-v0a
 511     if (id==9) {id1=5; id2=9;}               // t0ac-v0ac
 512     if (id==10) {id1=11; id2=12;}            // v0c0-v0a0
 513     if (id==11) {id1=15; id2=16;}            // v0c1-v0a1
 514     if (id==12) {id1=19; id2=20;}            // v0c2-v0a2
 515     if (id==13) {id1=23; id2=24;}            // v0c3-v0a3
 516     if (id==14) {id1=3; id2=12;}             // t0c-v0a0
 517     if (id==15) {id1=3; id2=16;}             // t0c-v0a1
 518     if (id==16) {id1=3; id2=20;}             // t0c-v0a2
 519     if (id==17) {id1=3; id2=24;}             // t0c-v0a3
 520     if (id==18) {id1=4; id2=11;}             // t0a-v0c0
 521     if (id==19) {id1=4; id2=15;}             // t0a-v0c1
 522     if (id==20) {id1=4; id2=19;}             // t0a-v0c2
 523     if (id==21) {id1=4; id2=23;}             // t0a-v0c3
 524     if (sumxy[ih1][id1][3]>-9000 && sumxy[ih2][id2][3]>-9000) {
 525      float dph=sumxy[ih1][id1][3]-sumxy[ih2][id2][3];
 526      res[ih][id][0]->Fill((float)icen,cos((ih+1.0)*dph));
 527      res[ih][id][1]->Fill((float)icen,sin((ih+1.0)*dph));
 528      float phi1=sumxy[ih1][id1][3];
 529      float phi2=sumxy[ih2][id2][3];
 530      if (ih>ih1) phi1*=(ih+1.0);
 531      else        phi1*=(ih1+1.0);
 532      if (ih>ih2) phi2*=(ih+1.0);
 533      else        phi2*=(ih2+1.0);
 534      phi1=atan2(sin(phi1),cos(phi1));
 535      phi2=atan2(sin(phi2),cos(phi2));
 536      cor[ih][id][icen/2]->Fill(phi1,phi2);
 537     }
 538    }
 539   }
 540
 541   corr00->Fill(1.0*ntr,v1s+v2s);
 542   corr01->Fill(1.0*ntr,t1s+t2s);
 543   corr02->Fill(1.0*ntr,z1s+z2s);
 544   corr03->Fill(v1s+v2s,t1s+t2s);
 545   corr04->Fill(v1s+v2s,z1s+z2s);
 546   corr05->Fill(t1s+t2s,z1s+z2s);
 547   corr10->Fill(1.0*ntr,v0sm);
 548   corr11->Fill(1.0*ntr,t0sm);
 549   corr12->Fill(1.0*ntr,z0sm);
 550   corr13->Fill(v0sm,t0sm);
 551   corr14->Fill(v0sm,z0sm);
 552   corr15->Fill(t0sm,z0sm);
 553   if (ievent%1000==0) cout << "eve " << jentry << " " << ievent << " : "
 554                       << jev << " " << iev << " : " << ntr << " : "
 555                       << z1s+z2s << " " << t1s+t2s << " " << v1s+v2s << " : "
 556                       << z0sm << " " << t0sm << " " << v0sm << endl;
 557   ievent++;
 558  }
 559  }
 560  ofstream ofs;
 561  ofs.open("AnaEveMyTree.cal");
 562  for (int i=0; i<8; i++) {
 563   ofs << det0g->GetBinContent(i+1) << " ";
 564   if (i%8==7) ofs << endl;
 565   cout << det0g->GetBinContent(i+1) << " ";
 566   if (i%8==7) cout << endl;
 567  }
 568  for (int i=0; i<24; i++) {
 569   ofs << det1g->GetBinContent(i+1) << " ";
 570   if (i%8==7) ofs << endl;
 571   cout << det1g->GetBinContent(i+1) << " ";
 572   if (i%8==7) cout << endl;
 573  }
 574  for (int i=0; i<64; i++) {
 575   ofs << det2g->GetBinContent(i+1) << " ";
 576   if (i%8==7) ofs << endl;
 577   cout << det2g->GetBinContent(i+1) << " ";
 578   if (i%8==7) cout << endl;
 579  }
 580  for (int ic=0; ic<10; ic++) {
 581   for (int iz=0; iz<10; iz++) {
 582    for (int ih=0; ih<4; ih++) {
 583     for (int id=0; id<27; id++) {
 584      for (int ib=0; ib<2; ib++) {
 585       ofs << ave[ic][iz][ih][id]->GetBinContent(ib+1) << " ";
 586       ofs << ave[ic][iz][ih][id]->GetBinError  (ib+1) << " ";
 587       if (ic==4 && iz==4 && ih==0) {
 588        cout << ave[ic][iz][ih][id]->GetBinContent(ib+1) << " ";
 589        cout << ave[ic][iz][ih][id]->GetBinError  (ib+1) << " ";
 590       }
 591      }
 592      ofs << endl;
 593      if (ic==4 && iz==4 && ih==0) cout << endl;
 594      for (int ib=0; ib<2; ib++) {
 595       for (int io=0; io<8; io++) {
 596        ofs << flt[ic][iz][ih][id]->GetBinContent(ib*8+io+1) << " ";
 597        if (ic==4 && iz==4 && ih==0) {
 598         cout << flt[ic][iz][ih][id]->GetBinContent(ib*8+io+1) << " ";
 599        }
 600       }
 601       ofs << endl;
 602       if (ic==4 && iz==4 && ih==0) cout << endl;
 603      }
 604     }
 605    }
 606   }
 607  }
 608  ofs.close();
 609  tf->Write();
 610  tf->Close();
 611 }