StRoot: StHbtMaker/CorrFctn/ParityDevCorrFctn.cxx Source File

00001 #include "StHbtMaker/CorrFctn/ParityDevCorrFctn.h"
00002 #include "StHbtMaker/Infrastructure/StHbtAnalysis.h"
00003 #include "StHbtMaker/Cut/mikesEventCut.h"
00004 
00005 #include <cstdio>
00006 
00007 #ifdef __ROOT__
00008 ClassImp(ParityDevCorrFctn) 
00009 #endif
00010 
00011 //____________________________
00012 ParityDevCorrFctn::ParityDevCorrFctn(){
00013   //mTagWriter = StHbtTagWriter::Instance();  // get the singleton
00014 
00015   // histograms
00016   mKMiSamex = new StHbt1DHisto("SameKMinusX","Parity Same KMinus x",100,-0.4,0.4);
00017   mKMiSamey = new StHbt1DHisto("SameKMinusY","Parity Same KMinus y",100,-0.4,0.4);
00018   mKMiSamez = new StHbt1DHisto("SameKMinusZ","Parity Same KMinus z",100,-0.4,0.4);
00019   mJcSamex = new StHbt1DHisto("SameJcX","Parity Same Jc X",100,-0.2,0.2);
00020   mJcSamey = new StHbt1DHisto("SameJcY","Parity Same Jc Y",100,-0.2,0.2);
00021   mJcSamez = new StHbt1DHisto("SameJcZ","Parity Same Jc Z",100,-0.2,0.2);
00022   mJcKtSame = new StHbt1DHisto("SameJcKt","Parity Same JcKt",100,-0.005,0.005);
00023   mKtwistKtSame = new StHbt1DHisto("SameKtwistKt","Parity Same KtwistKt",100,-0.005,0.005);
00024 
00025   mKMiMixedx = new StHbt1DHisto("MixedKMinusX","Parity Mixed KMinus x",100,-0.4,0.4);
00026   mKMiMixedy = new StHbt1DHisto("MixedKMinusY","Parity Mixed KMinus y",100,-0.4,0.4);
00027   mKMiMixedz = new StHbt1DHisto("MixedKMinusZ","Parity Mixed KMinus z",100,-0.4,0.4);
00028   mJcMixedx = new StHbt1DHisto("MixedJcX","Parity Mixed Jc X",100,-0.2,0.2);
00029   mJcMixedy = new StHbt1DHisto("MixedJcY","Parity Mixed Jc Y",100,-0.2,0.2);
00030   mJcMixedz = new StHbt1DHisto("MixedJcZ","Parity Mixed Jc Z",100,-0.2,0.2);
00031   mJcKtMixed = new StHbt1DHisto("MixedJcKt","Parity Mixed JcKt",100,-0.005,0.005);
00032   mKtwistKtMixed = new StHbt1DHisto("MixedKtwistKt","Parity Mixed KtwistKt",100,-0.005,0.005);
00033 
00034   mJcKtBinomial = new StHbt1DHisto("JcKtBin","Number JcKt positive",100,-0.10,0.10);
00035   mNumPairsBin = new StHbt1DHisto("Numpairs","Number of Pairs",2000,000.,2000.);
00036 
00037 
00038   // end histograms
00039 
00040   // this next bit is unfortunately needed so that we can have many histos of same "title"
00041   // it is neccessary if we typedef StHbt1DHisto to TH1d (which we do)
00042   mKMiSamex->SetDirectory(0);
00043   mKMiSamey->SetDirectory(0);
00044   mKMiSamez->SetDirectory(0);
00045   mJcSamex->SetDirectory(0);
00046   mJcSamey->SetDirectory(0);
00047   mJcSamez->SetDirectory(0);
00048   mJcKtSame->SetDirectory(0);
00049   mKtwistKtSame->SetDirectory(0);
00050   mKMiMixedx->SetDirectory(0);
00051   mKMiMixedy->SetDirectory(0);
00052   mKMiMixedz->SetDirectory(0);
00053   mJcMixedx->SetDirectory(0);
00054   mJcMixedy->SetDirectory(0);
00055   mJcMixedz->SetDirectory(0);
00056   mJcKtMixed->SetDirectory(0);
00057   mKtwistKtMixed->SetDirectory(0);
00058   mJcKtBinomial->SetDirectory(0);
00059   mNumPairsBin->SetDirectory(0); 
00060 
00061   mKMiMixedx->Sumw2();
00062   mKMiMixedy->Sumw2();
00063   mKMiMixedz->Sumw2();
00064   mJcMixedx->Sumw2();
00065   mJcMixedy->Sumw2();
00066   mJcMixedz->Sumw2();
00067   mJcKtMixed->Sumw2();
00068   mKtwistKtMixed->Sumw2();
00069   mKMiSamex->Sumw2();
00070   mKMiSamey->Sumw2();
00071   mKMiSamez->Sumw2();
00072   mJcSamex->Sumw2();
00073   mJcSamey->Sumw2();
00074   mJcSamez->Sumw2();
00075   mJcKtSame->Sumw2();
00076   mKtwistKtSame->Sumw2();
00077   mJcKtBinomial->Sumw2();
00078   mNumPairsBin->Sumw2();
00079 
00080 }   
00081  
00082 //____________________________
00083 ParityDevCorrFctn::~ParityDevCorrFctn(){
00084   // histograms
00085   delete mKMiSamex;
00086   delete mKMiSamey;
00087   delete mKMiSamez;
00088   delete mJcSamex ;
00089   delete mJcSamey ;
00090   delete mJcSamez ;
00091   delete mJcKtSame ;
00092   delete mKtwistKtSame ;
00093 
00094   delete mKMiMixedx;
00095   delete mKMiMixedy;
00096   delete mKMiMixedz;
00097   delete mJcMixedx ;
00098   delete mJcMixedy ;
00099   delete mJcMixedz ;
00100   delete mJcKtMixed ;
00101   delete mKtwistKtMixed ;
00102   // end histograms;
00103 }
00104 //_________________________
00105 void ParityDevCorrFctn::Finish(){
00106   cout << " alive in finish " << endl;
00107 
00108 }    
00109 //____________________________
00110 StHbtString ParityDevCorrFctn::Report(){
00111   string stemp = "Parity Correlation Function Report:\n";
00112   StHbtString returnThis = stemp;
00113   return returnThis;
00114 }
00115 //____________________________
00116 inline void ParityDevCorrFctn::AddRealPair(const StHbtPair* pair){
00117   cout << "WARNING ParityDevCorrFctn::AddRealPair shouldn't be called" << endl;
00118 }
00119 //_________________________
00120 
00121 //_________________________
00122 inline void ParityDevCorrFctn::AddMixedPair(const StHbtPair* pair){
00123   cout << "WARNING ParityDevCorrFctn::AddMixedPair shouldn't be called" << endl;
00124 }
00125 
00126 //_________________________
00127 inline void ParityDevCorrFctn::ParityCompute(ParityBuff *Plus, ParityBuff *Minus, int mxd){
00128 
00129       StHbtLorentzVector    PlusTrack; 
00130       StHbtLorentzVector    MinusTrack; 
00131       StHbtThreeVector      KMinusPos;
00132       StHbtThreeVector      KMinusNeg;
00133       StHbtThreeVector      KMinus;
00134       StHbtThreeVector      Kt;
00135 
00136       double  JcKt =0;
00137       double  PairJcKt = 0;
00138       double  JcKtBinom = 0;
00139       double  KtwistKt = 0;      
00140 
00141     int plusSize = Plus->size();
00142     int minusSize = Minus->size();
00143 
00144      if (mxd == SAME) {
00145          cout << "********  we got to the ParityCompute for SAME event" << endl;
00146          cout << " size of same event plus is "  << Plus->size() << "  ";
00147          cout << " size of same event minus is " << Minus->size() << endl;
00148       }
00149       if (mxd == MIXED) {
00150          cout << "********  we got to the ParityCompute for MIXED event" << endl;
00151          cout << " size of mixed event plus is "  << Plus->size() << "  ";
00152          cout << " size of mixed event minus is " << Minus->size() << endl;
00153       }
00154     //     ****** Compute K- **********
00155 
00156        {for (int jjj = 0; jjj <plusSize; jjj++){
00157          KMinusPos += (*Plus)[jjj].vect().unit(); 
00158        }}
00159        {for (int jjj = 0; jjj < minusSize; jjj++){
00160          KMinusNeg += (*Minus)[jjj].vect().unit(); 
00161        }}
00162        KMinus = (KMinusPos/plusSize) - (KMinusNeg/minusSize);
00163        Kt = KMinus;
00164        Kt.setZ(0.);
00165        cout << "Kt is " << Kt.x() << "," << Kt.y() << endl ; 
00166 
00167        // now calculate Ktwist*Kt
00168 
00169        double plusSum = 0.;
00170        {for (int jjj = 0; jjj < plusSize; jjj++){
00171          StHbtThreeVector TempV = (*Plus)[jjj].vect();
00172          double step1 = (TempV.y() * Kt.x()-TempV.x() * Kt.y() );  // this is y-component rotated to Kt coordinates
00173          double step2 = step1*TempV.z();                                                                     // multiply by z-component
00174          double step3 = step2/(TempV.mag2() * Kt.mag());                         // normalize
00175          plusSum += step3;
00176        }}
00177        double minusSum = 0.;
00178        {for (int jjj = 0; jjj < minusSize; jjj++){
00179          StHbtThreeVector TempV = (*Minus)[jjj].vect();
00180          double step1 = (TempV.y() * Kt.x()-TempV.x() * Kt.y() );  // this is y-component rotated to Kt coordinates
00181          double step2 = step1*TempV.z();                                                                     // multiply by z-component
00182          double step3 = step2/(TempV.mag2() * Kt.mag());                         // normalize
00183          minusSum += step3;
00184        }}
00185 
00186        double Ktwist      = (plusSum/ plusSize) - (minusSum/ minusSize);
00187        KtwistKt = Ktwist*Kt.mag();       
00188 
00189        cout << "KtwistKt  is " << KtwistKt << endl;
00190        // end Ktwist
00191 
00192        // now JcKt -here we use  unique pairs since it's approx. as sensitive as using all pairs (?)
00193 
00194        int smallSize = plusSize;
00195        int numPairs = 0;
00196        int numJcKtPlus = 0;
00197        if (minusSize < plusSize) smallSize = minusSize;
00198        StHbtThreeVector Jc;
00199         {for (int jjj = 0; jjj < smallSize; jjj++){
00200          PlusTrack = (*Plus)[jjj];
00201          MinusTrack = (*Minus)[jjj];
00202               StHbtThreeVector vOne  = PlusTrack.vect();
00203               StHbtThreeVector vTwo  = MinusTrack.vect();
00204                          StHbtThreeVector PairJc = vOne.unit().cross(vTwo.unit());
00205                           if ( (vOne.z()*vTwo.z()) < 0.){
00206                            PairJc.setX(-PairJc.x());
00207                }
00208                           if ( (vOne.x()*vTwo.x()) < 0.){
00209                            PairJc.setY(-PairJc.y());
00210                           }
00211                           if ( (vOne.y()*vTwo.y()) < 0.){
00212                            PairJc.setZ(-PairJc.z());
00213                           }
00214                Jc += PairJc;
00215                //this is for binomial only
00216                numPairs++;
00217                PairJcKt = PairJc.dot(Kt);
00218                if (PairJcKt > 0){
00219                  numJcKtPlus++;
00220                }
00221                // end binomial
00222        }}
00223        Jc = Jc/(smallSize); // remember to change normalization if we go back to all pairs
00224        JcKt = Jc.dot(Kt);
00225        cout << "JcKt  is " << JcKt << endl;
00226        JcKtBinom = ( double(numJcKtPlus)/double(numPairs) ) -1./2. ; 
00227        cout << "JcKt bin = " << numJcKtPlus << "/" <<numPairs << "-1/2  ="<< JcKtBinom << endl; 
00228       // end JcKt
00229 
00230       if (mxd == SAME) {
00231                    mKMiSamex->Fill(KMinus.x());
00232                    mKMiSamey->Fill(KMinus.y());
00233                    mKMiSamez->Fill(KMinus.z());
00234                    mJcSamex->Fill(Jc.x());
00235                    mJcSamey->Fill(Jc.y());
00236                    mJcSamez->Fill(Jc.z());
00237                    mJcKtSame->Fill(JcKt);
00238                    mKtwistKtSame->Fill(KtwistKt);
00239 
00240                    mJcKtBinomial->Fill(JcKtBinom); // binomial
00241                    mNumPairsBin->Fill(numPairs); // binomial               
00242       }
00243 
00244       if (mxd == MIXED) {
00245                    mKMiMixedx->Fill(KMinus.x());
00246                    mKMiMixedy->Fill(KMinus.y());
00247                    mKMiMixedz->Fill(KMinus.z());
00248                    mJcMixedx->Fill(Jc.x());
00249                    mJcMixedy->Fill(Jc.y());
00250                    mJcMixedz->Fill(Jc.z());
00251                    mJcKtMixed->Fill(JcKt);
00252                    mKtwistKtMixed->Fill(KtwistKt);
00253       }
00254 
00255 
00256 }
00257 //_________________________