PWG2/FEMTOSCOPY/AliFemtoUser/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 "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   AliFemtoCorrFctn(),
  20   fShareNumerator(0),
  21   fShareDenominator(0),
  22   fQualityNumerator(0),
  23   fQualityDenominator(0),
  24   fTPCSepNumerator(0),
  25   fTPCSepDenominator(0)
  26 {
  27   // set up numerator
  28   //  title = "Num Qinv (MeV/c)";
  29   char tTitNum[100] = "NumShare";
  30   strcat(tTitNum,title);
  31   fShareNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,1.00001);
  32   // set up denominator
  33   //title = "Den Qinv (MeV/c)";
  34   char tTitDen[100] = "DenShare";
  35   strcat(tTitDen,title);
  36   fShareDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,1.00001);
  37
  38   char tTit2Num[100] = "NumQuality";
  39   strcat(tTit2Num,title);
  40   fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001);
  41   // set up denominator
  42   //title = "Den Qinv (MeV/c)";
  43   char tTit2Den[100] = "DenQuality";
  44   strcat(tTit2Den,title);
  45   fQualityDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001);
  46   // set up ratio
  47   //title = "Ratio Qinv (MeV/c)";
  48   // this next bit is unfortunately needed so that we can have many histos of same "title"
  49   // it is neccessary if we typedef TH2D to TH1d (which we do)
  50   //mShareNumerator->SetDirectory(0);
  51   //mShareDenominator->SetDirectory(0);
  52   //mRatio->SetDirectory(0);
  53
  54   char tTit3Num[100] = "NumTPCSep";
  55   strcat(tTit3Num,title);
  56   fTPCSepNumerator = new TH2D(tTit3Num,title,nbins,QinvLo,QinvHi,150,0.0,100.0);
  57   // set up denominator
  58   //title = "Den Qinv (MeV/c)";
  59   char tTit3Den[100] = "DenTPCSep";
  60   strcat(tTit3Den,title);
  61   fTPCSepDenominator = new TH2D(tTit3Den,title,nbins,QinvLo,QinvHi,150,0.0,100.0);
  62
  63   // to enable error bar calculation...
  64   fShareNumerator->Sumw2();
  65   fShareDenominator->Sumw2();
  66
  67   fQualityNumerator->Sumw2();
  68   fQualityDenominator->Sumw2();
  69
  70   fTPCSepNumerator->Sumw2();
  71   fTPCSepDenominator->Sumw2();
  72 }
  73
  74 //____________________________
  75 AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(const AliFemtoShareQualityCorrFctn& aCorrFctn) :
  76   AliFemtoCorrFctn(),
  77   fShareNumerator(0),
  78   fShareDenominator(0),
  79   fQualityNumerator(0),
  80   fQualityDenominator(0),
  81   fTPCSepNumerator(0),
  82   fTPCSepDenominator(0)
  83 {
  84   // copy constructor
  85   if (aCorrFctn.fShareNumerator)
  86     fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
  87   if (aCorrFctn.fShareDenominator)
  88     fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
  89   if (aCorrFctn.fQualityNumerator)
  90     fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
  91   if (aCorrFctn.fQualityDenominator)
  92     fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
  93   if (aCorrFctn.fTPCSepNumerator)
  94     fTPCSepNumerator = new TH2D(*aCorrFctn.fTPCSepNumerator);
  95   if (aCorrFctn.fTPCSepDenominator)
  96     fTPCSepDenominator = new TH2D(*aCorrFctn.fTPCSepDenominator);
  97 }
  98 //____________________________
  99 AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
 100   // destructor
 101   delete fShareNumerator;
 102   delete fShareDenominator;
 103   delete fQualityNumerator;
 104   delete fQualityDenominator;
 105   delete fTPCSepNumerator;
 106   delete fTPCSepDenominator;
 107 }
 108 //_________________________
 109 AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
 110 {
 111   // assignment operator
 112   if (this == &aCorrFctn)
 113     return *this;
 114
 115   if (aCorrFctn.fShareNumerator)
 116     fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
 117   else
 118     fShareNumerator = 0;
 119   if (aCorrFctn.fShareDenominator)
 120     fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
 121   else
 122     fShareDenominator = 0;
 123   if (aCorrFctn.fQualityNumerator)
 124     fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
 125   else
 126     fQualityNumerator = 0;
 127   if (aCorrFctn.fQualityDenominator)
 128     fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
 129   else
 130     fQualityDenominator = 0;
 131   if (aCorrFctn.fTPCSepNumerator)
 132     fTPCSepNumerator = new TH2D(*aCorrFctn.fTPCSepNumerator);
 133   else
 134     fTPCSepNumerator = 0;
 135   if (aCorrFctn.fTPCSepDenominator)
 136     fTPCSepDenominator = new TH2D(*aCorrFctn.fTPCSepDenominator);
 137   else
 138     fTPCSepDenominator = 0;
 139
 140   return *this;
 141 }
 142 //_________________________
 143 void AliFemtoShareQualityCorrFctn::Finish(){
 144   // here is where we should normalize, fit, etc...
 145   // we should NOT Draw() the histos (as I had done it below),
 146   // since we want to insulate ourselves from root at this level
 147   // of the code.  Do it instead at root command line with browser.
 148   //  mShareNumerator->Draw();
 149   //mShareDenominator->Draw();
 150   //mRatio->Draw();
 151
 152 }
 153
 154 //____________________________
 155 AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
 156   // create report
 157   string stemp = "Qinv Correlation Function Report:\n";
 158   char ctemp[100];
 159   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
 160   stemp += ctemp;
 161   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
 162   stemp += ctemp;
 163   //  stemp += mCoulombWeight->Report();
 164   AliFemtoString returnThis = stemp;
 165   return returnThis;
 166 }
 167 //____________________________
 168 void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){
 169   // add real (effect) pair
 170   if (fPairCut)
 171     if (!fPairCut->Pass(pair)) return;
 172
 173   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 174   Int_t nh = 0;
 175   Int_t an = 0;
 176   Int_t ns = 0;
 177
 178   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
 179     // If both have clusters in the same row
 180     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
 181         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 182       // Do they share it ?
 183       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
 184           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
 185         {
 186 //        if (tQinv < 0.01) {
 187 //          cout << "Shared cluster in row " << imap << endl;
 188 //        }
 189           an++;
 190           nh+=2;
 191           ns+=2;
 192         }
 193
 194       // Different hits on the same padrow
 195       else {
 196         an--;
 197         nh+=2;
 198       }
 199     }
 200     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
 201              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 202       // One track has a hit, the other does not
 203       an++;
 204       nh++;
 205     }
 206   }
 207 //    if (tQinv < 0.01) {
 208 //     cout << "Qinv of the pair is " << tQinv << endl;
 209 //     cout << "Clusters: " << endl;
 210 //     for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
 211 //       cout << imap ;
 212 //       if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
 213 //       else cout << " 0 " ;
 214 //       if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
 215 //       else cout << " 0 " ;
 216 //       cout << "     ";
 217 //       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
 218 //       else cout << " X ";
 219 //       if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
 220 //       else cout << " X ";
 221 //       cout << endl;
 222 //     }
 223 //   }
 224
 225   Float_t hsmval = 0.0;
 226   Float_t hsfval = 0.0;
 227
 228   if (nh >0) {
 229     hsmval = an*1.0/nh;
 230     hsfval = ns*1.0/nh;
 231   }
 232
 233 //   if ((tQinv < 0.005) && (hsmval<-0.0)) {
 234 //     cout << "Quality  Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl;
 235 //     cout << "Qinv of the pair is " << tQinv << endl;
 236 //     cout << "Clusters: " << endl;
 237 //     for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
 238 //       cout << imap ;
 239 //       if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
 240 //       else cout << " 0 " ;
 241 //       if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
 242 //       else cout << " 0 " ;
 243 //       cout << "     ";
 244 //       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
 245 //       else cout << " X ";
 246 //       if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
 247 //       else cout << " X ";
 248 //       cout << endl;
 249 //     }
 250 //     cout << "Momentum1 "
 251 //       << pair->Track1()->Track()->P().x() << " "
 252 //       << pair->Track1()->Track()->P().y() << " "
 253 //       << pair->Track1()->Track()->P().z() << " "
 254 //       << pair->Track1()->Track()->Label() << " "
 255 //       << pair->Track1()->Track()->TrackId() << " "
 256 //       << pair->Track1()->Track()->Flags() << " "
 257 //       << pair->Track1()->Track()->KinkIndex(0) << " "
 258 //       << pair->Track1()->Track()->KinkIndex(1) << " "
 259 //       << pair->Track1()->Track()->KinkIndex(2) << " "
 260 //       << pair->Track1()->Track()->ITSchi2() << " "
 261 //       << pair->Track1()->Track()->ITSncls() << " "
 262 //       << pair->Track1()->Track()->TPCchi2() << " "
 263 //       << pair->Track1()->Track()->TPCncls() << " "
 264 //       << endl;
 265 //     cout << "Momentum2 "
 266 //       << pair->Track2()->Track()->P().x() << " "
 267 //       << pair->Track2()->Track()->P().y() << " "
 268 //       << pair->Track2()->Track()->P().z() << " "
 269 //       << pair->Track2()->Track()->Label() << " "
 270 //       << pair->Track2()->Track()->TrackId() << " "
 271 //       << pair->Track2()->Track()->Flags() << " "
 272 //       << pair->Track2()->Track()->KinkIndex(0) << " "
 273 //       << pair->Track2()->Track()->KinkIndex(1) << " "
 274 //       << pair->Track2()->Track()->KinkIndex(2) << " "
 275 //       << pair->Track2()->Track()->ITSchi2() << " "
 276 //       << pair->Track2()->Track()->ITSncls() << " "
 277 //       << pair->Track2()->Track()->TPCchi2() << " "
 278 //       << pair->Track2()->Track()->TPCncls() << " "
 279 //       << endl;
 280 //   }
 281
 282   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 283   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 284   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 285   double dist = sqrt(distx*distx + disty*disty + distz*distz);
 286
 287   fShareNumerator->Fill(tQinv, hsfval);
 288   fQualityNumerator->Fill(tQinv, hsmval);
 289   fTPCSepNumerator->Fill(tQinv, dist);
 290
 291   //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
 292   //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
 293 }
 294 //____________________________
 295 void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
 296   // add mixed (background) pair
 297   if (fPairCut)
 298     if (!fPairCut->Pass(pair)) return;
 299
 300   double weight = 1.0;
 301   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 302   Int_t nh = 0;
 303   Int_t an = 0;
 304   Int_t ns = 0;
 305
 306   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
 307     // If both have clusters in the same row
 308     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
 309         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 310       // Do they share it ?
 311       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
 312           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
 313         {
 314           //      cout << "A shared cluster !!!" << endl;
 315           //    cout << "imap idx1 idx2 "
 316           //         << imap << " "
 317           //         << tP1idx[imap] << " " << tP2idx[imap] << endl;
 318           an++;
 319           nh+=2;
 320           ns+=2;
 321         }
 322
 323       // Different hits on the same padrow
 324       else {
 325         an--;
 326         nh+=2;
 327       }
 328     }
 329     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
 330              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
 331       // One track has a hit, the other does not
 332       an++;
 333       nh++;
 334     }
 335   }
 336
 337   Float_t hsmval = 0.0;
 338   Float_t hsfval = 0.0;
 339
 340   if (nh >0) {
 341     hsmval = an*1.0/nh;
 342     hsfval = ns*1.0/nh;
 343   }
 344
 345   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
 346   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
 347   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
 348   double dist = sqrt(distx*distx + disty*disty + distz*distz);
 349
 350   fShareDenominator->Fill(tQinv,hsfval,weight);
 351   fQualityDenominator->Fill(tQinv,hsmval,weight);
 352   fTPCSepDenominator->Fill(tQinv, dist, weight);
 353
 354 }
 355
 356
 357 void AliFemtoShareQualityCorrFctn::WriteHistos()
 358 {
 359   // Write out result histograms
 360   fShareNumerator->Write();
 361   fShareDenominator->Write();
 362   fQualityNumerator->Write();
 363   fQualityDenominator->Write();
 364   fTPCSepNumerator->Write();
 365   fTPCSepDenominator->Write();
 366
 367 }
 368 //______________________________
 369 TList* AliFemtoShareQualityCorrFctn::GetOutputList()
 370 {
 371   // Prepare the list of objects to be written to the output
 372   TList *tOutputList = new TList();
 373
 374   tOutputList->Add(fShareNumerator);
 375   tOutputList->Add(fShareDenominator);
 376   tOutputList->Add(fQualityNumerator);
 377   tOutputList->Add(fQualityDenominator);
 378   tOutputList->Add(fTPCSepNumerator);
 379   tOutputList->Add(fTPCSepDenominator);
 380
 381   return tOutputList;
 382 }