PWG2/FEMTOSCOPY/AliFemtoUser/CorrFctn/AliFemtoShareQualityCorrFctn.cxx

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