PWG2/FEMTOSCOPY/AliFemtoUser/CorrFctn/AliFemtoShareQualityCorrFctn.cxx

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
   6  ***************************************************************************
   7  *
   8  * Description: part of STAR HBT Framework: AliFemtoMaker package
   9  *   a simple Q-invariant correlation function
  10  *
  11  ***************************************************************************
  12  *
  13  * $Log$
  14  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  15  * First version on CVS
  16  *
  17  * Revision 1.4  2000/01/25 17:34:45  laue
  18  * I. In order to run the stand alone version of the AliFemtoMaker the following
  19  * changes have been done:
  20  * a) all ClassDefs and ClassImps have been put into #ifdef __ROOT__ statements
  21  * b) unnecessary includes of StMaker.h have been removed
  22  * c) the subdirectory AliFemtoMaker/doc/Make has been created including everything
  23  * needed for the stand alone version
  24  *
  25  * II. To reduce the amount of compiler warning
  26  * a) some variables have been type casted
  27  * b) some destructors have been declared as virtual
  28  *
  29  * Revision 1.3  1999/07/29 02:47:09  lisa
  30  * 1) add OpeningAngle correlation function 2) add AliFemtoMcEventReader 3) make histos in CorrFctns do errors correctly
  31  *
  32  * Revision 1.2  1999/07/06 22:33:20  lisa
  33  * Adjusted all to work in pro and new - dev itself is broken
  34  *
  35  * Revision 1.1.1.1  1999/06/29 16:02:57  lisa
  36  * Installation of AliFemtoMaker
  37  *
  38  **************************************************************************/
  39
  40 #include "CorrFctn/AliFemtoShareQualityCorrFctn.h"
  41 //#include "Infrastructure/AliFemtoHisto.hh"
  42 #include <cstdio>
  43
  44 #ifdef __ROOT__
  45 ClassImp(AliFemtoShareQualityCorrFctn)
  46 #endif
  47
  48 //____________________________
  49 AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi){
  50   // set up numerator
  51   //  title = "Num Qinv (MeV/c)";
  52   char TitNum[100] = "NumShare";
  53   strcat(TitNum,title);
  54   fShareNumerator = new TH2D(TitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
  55   // set up denominator
  56   //title = "Den Qinv (MeV/c)";
  57   char TitDen[100] = "DenShare";
  58   strcat(TitDen,title);
  59   fShareDenominator = new TH2D(TitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
  60
  61   char Tit2Num[100] = "NumQuality";
  62   strcat(Tit2Num,title);
  63   fQualityNumerator = new TH2D(Tit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
  64   // set up denominator
  65   //title = "Den Qinv (MeV/c)";
  66   char Tit2Den[100] = "DenQuality";
  67   strcat(Tit2Den,title);
  68   fQualityDenominator = new TH2D(Tit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
  69   // set up ratio
  70   //title = "Ratio Qinv (MeV/c)";
  71   // this next bit is unfortunately needed so that we can have many histos of same "title"
  72   // it is neccessary if we typedef TH2D to TH1d (which we do)
  73   //mShareNumerator->SetDirectory(0);
  74   //mShareDenominator->SetDirectory(0);
  75   //mRatio->SetDirectory(0);
  76
  77   // to enable error bar calculation...
  78   fShareNumerator->Sumw2();
  79   fShareDenominator->Sumw2();
  80
  81   fQualityNumerator->Sumw2();
  82   fQualityDenominator->Sumw2();
  83 }
  84
  85 //____________________________
  86 AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
  87   delete fShareNumerator;
  88   delete fShareDenominator;
  89   delete fQualityNumerator;
  90   delete fQualityDenominator;
  91 }
  92 //_________________________
  93 void AliFemtoShareQualityCorrFctn::Finish(){
  94   // here is where we should normalize, fit, etc...
  95   // we should NOT Draw() the histos (as I had done it below),
  96   // since we want to insulate ourselves from root at this level
  97   // of the code.  Do it instead at root command line with browser.
  98   //  mShareNumerator->Draw();
  99   //mShareDenominator->Draw();
 100   //mRatio->Draw();
 101
 102 }
 103
 104 //____________________________
 105 AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
 106   string stemp = "Qinv Correlation Function Report:\n";
 107   char ctemp[100];
 108   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
 109   stemp += ctemp;
 110   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
 111   stemp += ctemp;
 112   //  stemp += mCoulombWeight->Report();
 113   AliFemtoString returnThis = stemp;
 114   return returnThis;
 115 }
 116 //____________________________
 117 void AliFemtoShareQualityCorrFctn::AddRealPair(const AliFemtoPair* pair){
 118   double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
 119   Int_t nh = 0;
 120   Int_t an = 0;
 121   Int_t ns = 0;
 122
 123   for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
 124     // If both have clusters in the same row
 125     if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
 126         pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 127       // Do they share it ?
 128       if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) ||
 129           pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
 130         {
 131           if (Qinv < 0.01) {
 132             cout << "Shared cluster in row " << imap << endl;
 133           }
 134           an++;
 135           nh+=2;
 136           ns+=2;
 137         }
 138
 139       // Different hits on the same padrow
 140       else {
 141         an--;
 142         nh+=2;
 143       }
 144     }
 145     else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) ||
 146              pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 147       // One track has a hit, the other does not
 148       an++;
 149       nh++;
 150     }
 151   }
 152   if (Qinv < 0.01) {
 153     cout << "Qinv of the pair is " << Qinv << endl;
 154     cout << "Clusters: " << endl;
 155     for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
 156       cout << imap ;
 157       if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
 158       else cout << " 0 " ;
 159       if (pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
 160       else cout << " 0 " ;
 161       cout << "     ";
 162       if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
 163       else cout << " X ";
 164       if (pair->track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
 165       else cout << " X ";
 166       cout << endl;
 167     }
 168   }
 169
 170   Float_t hsmval = 0.0;
 171   Float_t hsfval = 0.0;
 172
 173   if (nh >0) {
 174     hsmval = an*1.0/nh;
 175     hsfval = ns*1.0/nh;
 176   }
 177
 178   if (Qinv < 0.01) {
 179     cout << "Quality  Sharity " << hsmval << " " << hsfval << " " << pair->track1()->Track() << " " << pair->track2()->Track() << endl;
 180   }
 181
 182   fShareNumerator->Fill(Qinv, hsfval);
 183   fQualityNumerator->Fill(Qinv, hsmval);
 184   //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
 185   //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
 186 }
 187 //____________________________
 188 void AliFemtoShareQualityCorrFctn::AddMixedPair(const AliFemtoPair* pair){
 189   double weight = 1.0;
 190   double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
 191   Int_t nh = 0;
 192   Int_t an = 0;
 193   Int_t ns = 0;
 194
 195   for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
 196     // If both have clusters in the same row
 197     if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
 198         pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 199       // Do they share it ?
 200       if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) ||
 201           pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
 202         {
 203           //      cout << "A shared cluster !!!" << endl;
 204           //    cout << "imap idx1 idx2 "
 205           //         << imap << " "
 206           //         << tP1idx[imap] << " " << tP2idx[imap] << endl;
 207           an++;
 208           nh+=2;
 209           ns+=2;
 210         }
 211
 212       // Different hits on the same padrow
 213       else {
 214         an--;
 215         nh+=2;
 216       }
 217     }
 218     else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) ||
 219              pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 220       // One track has a hit, the other does not
 221       an++;
 222       nh++;
 223     }
 224   }
 225
 226   Float_t hsmval = 0.0;
 227   Float_t hsfval = 0.0;
 228
 229   if (nh >0) {
 230     hsmval = an*1.0/nh;
 231     hsfval = ns*1.0/nh;
 232   }
 233
 234   fShareDenominator->Fill(Qinv,hsfval,weight);
 235   fQualityDenominator->Fill(Qinv,hsmval,weight);
 236 }
 237
 238
 239 void AliFemtoShareQualityCorrFctn::WriteHistos()
 240 {
 241   fShareNumerator->Write();
 242   fShareDenominator->Write();
 243   fQualityNumerator->Write();
 244   fQualityDenominator->Write();
 245
 246 }