PWGCF/FEMTOSCOPY/AliFemto/AliTwoTrackRes.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //345678901234567890123456789012345678901234567890123456789012345678901234567890
   3 //       1         2         3         4         5         6         7         8
   4 //
   5 // Tool to study two-track effects in ALICE for femtoscopic analyses
   6 // J. Mercado <mercado@physi.uni-heidelberg.de> Last modified: 20.01.2011
   7 //
   8 ////////////////////////////////////////////////////////////////////////////////
   9
  10 #include "AliTwoTrackRes.h"
  11 #include <iostream>
  12 #include <fstream>
  13 #include "TMath.h"
  14 #include "TROOT.h"
  15 #include "TFile.h"
  16 #include "TChain.h"
  17 #include "TLeaf.h"
  18 #include "TNtuple.h"
  19 #include "TRandom2.h"
  20
  21 ClassImp(AliTwoTrackRes)
  22
  23 //______________________________________________________________________________
  24 // Constructor(s)
  25
  26 AliTwoTrackRes::AliTwoTrackRes(const char *name) :
  27   AliAnalysisTask(name,""), fChain(0), fESDEvent(0),
  28   fOutContainer(0), fTrackCuts(0), fNTuple1(0),
  29   fNTuple2(0), fP1(), fP2(), fPb1(), fPb2(), fP(), fQ(), fTpcEnt1(), fTpcEnt2(),
  30   fTpcDist(), fOutFilename()
  31 {
  32   DefineInput(0, TChain::Class());     // Slot input 0 reads from a TChain
  33   DefineOutput(0, TObjArray::Class()); // Slot output 0 writes into a TObjArray
  34 }
  35
  36 AliTwoTrackRes::AliTwoTrackRes(const AliTwoTrackRes& aTwoTrackRes) :
  37   AliAnalysisTask(aTwoTrackRes), fChain(0), fESDEvent(0), fOutContainer(0),
  38   fTrackCuts(0), fNTuple1(0), fNTuple2(0), fP1(), fP2(), fPb1(), fPb2(), fP(),
  39   fQ(), fTpcEnt1(), fTpcEnt2(), fTpcDist(), fOutFilename()
  40 {
  41   //Copy constructor
  42   fChain = aTwoTrackRes.fChain;
  43   fESDEvent = aTwoTrackRes.fESDEvent;
  44   fOutContainer = aTwoTrackRes.fOutContainer;
  45   fTrackCuts = aTwoTrackRes.fTrackCuts;
  46   fNTuple1 = aTwoTrackRes.fNTuple1;
  47   fNTuple2 = aTwoTrackRes.fNTuple2;
  48   fP1 = aTwoTrackRes.fP1;
  49   fP2 = aTwoTrackRes.fP2;
  50   fPb1 = aTwoTrackRes.fPb1;
  51   fPb2 = aTwoTrackRes.fPb2;
  52   fP = aTwoTrackRes.fP;
  53   fQ = aTwoTrackRes.fQ;
  54   fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
  55   fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
  56   fTpcDist = aTwoTrackRes.fTpcDist;
  57   fOutFilename = aTwoTrackRes.fOutFilename;
  58 }
  59
  60 AliTwoTrackRes& AliTwoTrackRes::operator=(const AliTwoTrackRes& aTwoTrackRes)
  61 {
  62   // Assignment operator
  63   if (this == &aTwoTrackRes)
  64     return *this;
  65   fChain = aTwoTrackRes.fChain;
  66   fESDEvent = aTwoTrackRes.fESDEvent;
  67   fOutContainer = aTwoTrackRes.fOutContainer;
  68   fTrackCuts = aTwoTrackRes.fTrackCuts;
  69   fNTuple1 = aTwoTrackRes.fNTuple1;
  70   fNTuple2 = aTwoTrackRes.fNTuple2;
  71   fP1 = aTwoTrackRes.fP1;
  72   fP2 = aTwoTrackRes.fP2;
  73   fPb1 = aTwoTrackRes.fPb1;
  74   fPb2 = aTwoTrackRes.fPb2;
  75   fP = aTwoTrackRes.fP;
  76   fQ = aTwoTrackRes.fQ;
  77   fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
  78   fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
  79   fTpcDist = aTwoTrackRes.fTpcDist;
  80   fOutFilename = aTwoTrackRes.fOutFilename;
  81   return *this;
  82 }
  83
  84 AliTwoTrackRes::~AliTwoTrackRes() {printf("AliTwoTrackRes destroyed\n");}
  85
  86 void AliTwoTrackRes::ConnectInputData(Option_t *) {
  87 //______________________________________________________________________________
  88 // Connect input data and initialize track cuts
  89
  90   fChain = (TChain*)GetInputData(0);
  91   fESDEvent = new AliESDEvent();
  92   fESDEvent->ReadFromTree(fChain);
  93
  94   // Cuts to select primary tracks (ITS+TPC)
  95   fTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
  96   Double_t cov1, cov2, cov3, cov4, cov5; // diagonal cov. matrix elements
  97   Double_t nSigma;                       // max. DCA to primary vertex
  98   Int_t minNClustersTPC;                 // min. number of clusters per TPC tracks
  99   Double_t maxChi2PerClusterTPC;         // max. chi2 per cluster per TPC track
 100   Int_t cutMode = 1;                     // select cut mode
 101   if (cutMode == 1) {
 102   cov1 = 2; cov2 = 2; cov3 = 0.5; cov4 = 0.5; cov5 = 2;
 103   nSigma = 3; minNClustersTPC = 75; maxChi2PerClusterTPC = 3.5;
 104   fTrackCuts->SetMaxCovDiagonalElements(cov1, cov2, cov3, cov4, cov5);
 105   fTrackCuts->SetMaxNsigmaToVertex(nSigma);
 106   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 107   fTrackCuts->SetRequireTPCRefit(kTRUE);
 108   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 109   fTrackCuts->SetMinNClustersTPC(minNClustersTPC);
 110   fTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
 111   TString tag("Global tracking");}
 112 }
 113
 114 void AliTwoTrackRes::CreateOutputObjects() {
 115 //______________________________________________________________________________
 116 // Create output objects
 117
 118   fNTuple1 = new TNtuple("nt1","True pairs",
 119   "pt1:eta1:phi1:nsh1:pt2:eta2:phi2:nsh2:qinv:mindist:dist:corr:qfac");
 120   fNTuple2 = new TNtuple("nt2","Mixed pairs",
 121   "pt1:eta1:phi1:nsh1:pt2:eta2:phi2:nsh2:qinv:mindist:dist:corr:qfac");
 122   Int_t c = 0;
 123   fOutContainer = new TObjArray(2);
 124   fOutContainer->AddAt(fNTuple1, c++);
 125   fOutContainer->AddAt(fNTuple2, c++);
 126 }
 127
 128 void AliTwoTrackRes::Exec(Option_t *) {
 129 //______________________________________________________________________________
 130 // Create true and mixed pairs keeping some track parameters
 131
 132   double bfield = 5.0;
 133   static int nr=0;
 134   if (nr == 0) printf("\tStarting event loop...\n");
 135   printf("\rProcessing event %8d", nr);
 136   Double_t mpi = 0.13957; // [GeV/c^2]
 137   Double_t pidTrk1[AliPID::kSPECIES], pidTrk2[AliPID::kSPECIES];
 138   Int_t tpcIn = 80;   // [cm]
 139   Int_t tpcOut = 250; // [cm]
 140   Double_t tdist[170];
 141   Double_t tdistrot[170];
 142   Double_t tpcEnt1[3], tpcEnt2[3], pos1[3];
 143   TVector3 x1, x2, diff;
 144   TBits clu1, clu2, sha1, sha2;
 145   TRandom2 rnd;
 146   Int_t  ntracks = fESDEvent->GetNumberOfTracks();
 147   for(Int_t itrack = 0; itrack < ntracks; itrack++) {
 148     AliESDtrack *track1 = fESDEvent->GetTrack(itrack);
 149     AliExternalTrackParam *trp1 = const_cast<AliExternalTrackParam*>
 150       (track1->GetTPCInnerParam());
 151     if (!trp1) continue;
 152     if (!track1->IsOn(AliESDtrack::kTPCpid)) continue;
 153     track1->GetTPCpid(pidTrk1);
 154     Int_t q1 = trp1->Charge();
 155     if (!((fTrackCuts->AcceptTrack(track1)) && (q1 == 1) &&
 156           (pidTrk1[AliPID::kPion]+pidTrk1[AliPID::kMuon] > 0.5))) continue;
 157     if (!track1->GetInnerXYZ(tpcEnt1)) continue;
 158     clu1 = track1->GetTPCClusterMap();
 159     sha1 = track1->GetTPCSharedMap();
 160     SetTr1(track1->Pt(), track1->Eta(), track1->Phi(), mpi);
 161     SetTpcEnt1(tpcEnt1[0], tpcEnt1[1], tpcEnt1[2]);
 162     for(Int_t jtrack = 0; jtrack < itrack; jtrack++) {
 163       AliESDtrack *track2 = fESDEvent->GetTrack(jtrack);
 164       AliExternalTrackParam *trp2 = const_cast<AliExternalTrackParam*>
 165         (track2->GetTPCInnerParam());
 166       if (!trp2) continue;
 167       if (!track2->IsOn(AliESDtrack::kTPCpid)) continue;
 168       track2->GetTPCpid(pidTrk2);
 169       Int_t q2 = trp2->Charge();
 170       if (!((fTrackCuts->AcceptTrack(track2)) && (q2 == 1) &&
 171             (pidTrk2[AliPID::kPion]+pidTrk2[AliPID::kMuon] > 0.5))) continue;
 172       if (!track2->GetInnerXYZ(tpcEnt2)) continue;
 173       clu2 = track2->GetTPCClusterMap();
 174       sha2 = track2->GetTPCSharedMap();
 175       SetTr2(track2->Pt(), track2->Eta(), track2->Phi(), mpi);
 176       SetTpcEnt2(tpcEnt2[0], tpcEnt2[1], tpcEnt2[2]);
 177       for (Int_t i = tpcIn; i < tpcOut; i++) { // Minimum distance
 178         trp1->GetDistance(trp2, (double) i, pos1, bfield);
 179         x1.SetXYZ(pos1[0], pos1[1], pos1[2]);
 180         tdist[i-tpcIn] = x1.Mag();
 181         x1.SetXYZ(-pos1[0], -pos1[1], pos1[2]);
 182         tdistrot[i-tpcIn] = x1.Mag();
 183       }
 184       Double_t mindist = 100000;
 185       Int_t jmin=0;
 186       for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
 187         if (tdist[j] < mindist) {jmin=j;  mindist = tdist[j]; }
 188       }
 189       //      Double_t mindist = MinDist(track1, track2);
 190       Double_t dist = Dist();
 191       Double_t dphi = DPhi();
 192       Double_t deta = DEta();
 193       Int_t    nsh1 = GetNSha(clu1, sha1);
 194       Int_t    nsh2 = GetNSha(clu2, sha2);
 195       Double_t corr = Corr(clu1, clu2, sha1, sha2);
 196       Double_t qfac = Qfac(clu1, clu2, sha1, sha2);
 197       if ((TMath::Abs(track1->Eta())>0.8)&&(TMath::Abs(track2->Eta())>0.8)) continue;
 198       if ((TMath::Abs(dphi)<0.35)&&(deta<0.35)) {
 199       FillNTuple1(mindist,dist,corr,qfac,nsh1,nsh2);}    // True
 200       Double_t tr2rot = RotTr2Phi();                // Rotate trck2
 201       SetTr2(track2->Pt(), track2->Eta(), tr2rot, mpi);
 202       tpcEnt2[0] = -tpcEnt2[0];
 203       tpcEnt2[1] = -tpcEnt2[1];
 204       Double_t distrot = Dist();
 205       Double_t dphirot = DPhi();
 206       Double_t mindistrot = 100000;
 207       jmin=0;
 208       for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
 209         if (tdistrot[j] < mindistrot) {jmin=j;  mindistrot = tdistrot[j]; }
 210       }
 211       if ((TMath::Abs(dphirot)<0.35)&&(deta<0.35)) {
 212         if (rnd.Rndm() < 0.5) NoSwap();
 213         else Swap();
 214         FillNTuple2(mindistrot,distrot,corr,qfac,nsh1,nsh2);} // Mixed
 215     }
 216   }
 217   PostData(0, fOutContainer);
 218   nr++;
 219 }
 220
 221 void AliTwoTrackRes::Terminate(Option_t *) {
 222 //______________________________________________________________________________
 223 // Write output and clean up
 224
 225   fOutContainer = (TObjArray*)GetOutputData(0);
 226   TFile *f1  = new TFile( fOutFilename, "RECREATE" );
 227   fOutContainer->Write();
 228   f1->Flush();
 229   f1->Close();
 230   delete f1;
 231   delete fChain;
 232   delete fNTuple1;
 233   delete fNTuple2;
 234   printf("\n");
 235 }
 236
 237 //______________________________________________________________________________
 238 // Miscellaneous methods
 239
 240 // Set tracks
 241 void AliTwoTrackRes::SetTr1(double pt1, double eta1, double phi1, double m) {
 242   fP1.SetPtEtaPhiM(pt1, eta1, phi1, m);}
 243 void AliTwoTrackRes::SetTr2(double pt2, double eta2, double phi2, double m) {
 244   fP2.SetPtEtaPhiM(pt2, eta2, phi2, m);}
 245
 246 // Set nominal TPC entrance coordinates
 247 void AliTwoTrackRes::SetTpcEnt1(double x1, double y1, double z1) {
 248   fTpcEnt1.SetX(x1); fTpcEnt1.SetY(y1); fTpcEnt1.SetZ(z1);}
 249 void AliTwoTrackRes::SetTpcEnt2(double x2, double y2, double z2) {
 250   fTpcEnt2.SetX(x2); fTpcEnt2.SetY(y2); fTpcEnt2.SetZ(z2);}
 251
 252 double AliTwoTrackRes::MinDist(AliExternalTrackParam *trk1,
 253                                AliExternalTrackParam *trk2) {
 254 // Calculate minimum track separation within the TPC
 255
 256   int tpcIn = 0;   // [cm]
 257   int tpcOut = 170; // [cm]
 258   double tdist[170], pos[3];
 259   TVector3 x;
 260   for (int i = tpcIn; i < tpcOut; i++) {
 261     trk1->GetDistance(trk2, i, pos, 5000);
 262     x.SetXYZ(pos[0], pos[1], pos[2]);
 263     tdist[i-tpcIn] = x.Mag();
 264   }
 265   double maxdist = 0.0;
 266   for (int j = 0; j < tpcOut-tpcIn; j++) {
 267     if (tdist[j] > maxdist) { maxdist = tdist[j]; }
 268   }
 269   double mindist = maxdist;
 270   for (int j = 0; j < tpcOut-tpcIn; j++) {
 271     if (tdist[j] < mindist) { mindist = tdist[j]; }
 272   }
 273   return mindist;}
 274
 275 int AliTwoTrackRes::GetNSha(TBits cl, TBits sh) {
 276 // Get number of shared clusters
 277
 278   int ncl = cl.GetNbits();
 279   int sum = 0;
 280   for(int i = 0; i < ncl; i++) {
 281     if (!cl.TestBitNumber(i)) continue;
 282     int n = sh.TestBitNumber(i);
 283     sum += n;}
 284   return sum;}
 285
 286 double AliTwoTrackRes::Corr(TBits cl1,  TBits cl2, TBits sh1, TBits sh2) {
 287 // Calculate correlation coefficient
 288
 289   int ncl1 = cl1.GetNbits();
 290   int ncl2 = cl2.GetNbits();
 291   double sumN = 0;  double sumX = 0;  double sumY = 0;
 292   double sumXX = 0; double sumYY = 0; double sumXY = 0; double corr = -2.0;
 293   for(int i = 0; i < ncl1 && i < ncl2; i++) {
 294     if (!(cl1.TestBitNumber(i)&&cl2.TestBitNumber(i))) continue;
 295     int x = sh1.TestBitNumber(i);
 296     int y = sh2.TestBitNumber(i);
 297     sumN += 1.0;
 298     sumX += x;
 299     sumY += y;
 300     sumXX += x*x;
 301     sumYY += y*y;
 302     sumXY += x*y;
 303   }
 304   double meanX = sumX/sumN;
 305   double meanY = sumY/sumN;
 306   double meanXX = sumXX/sumN;
 307   double meanYY = sumYY/sumN;
 308   double meanXY = sumXY/sumN;
 309   double sX = TMath::Sqrt(TMath::Abs(meanXX-meanX*meanX));
 310   double sY = TMath::Sqrt(TMath::Abs(meanYY-meanY*meanY));
 311   if (sX*sY!=0) corr = (meanXY-meanX*meanY)/(sX*sY);
 312   return corr;}
 313
 314 double AliTwoTrackRes::Qfac(TBits cl1,  TBits cl2, TBits sh1, TBits sh2) {
 315 // Quality factor from AliFemto
 316
 317   int ncl1 = cl1.GetNbits();
 318   int ncl2 = cl2.GetNbits();
 319   int sumCls = 0; int sumSha = 0; int sumQ = 0;
 320   double shfrac = 0; double qfactor = 0;
 321   for(int i = 0; i < ncl1 && i < ncl2; i++) {
 322     if (cl1.TestBitNumber(i) && cl2.TestBitNumber(i)) { // Both clusters
 323       if (sh1.TestBitNumber(i) && sh2.TestBitNumber(i)) { // Shared
 324         sumQ++;
 325         sumCls+=2;
 326         sumSha+=2;}
 327       else {sumQ--; sumCls+=2;}
 328     }
 329     else if (cl1.TestBitNumber(i) || cl2.TestBitNumber(i)) { // Non shared
 330       sumQ++;
 331       sumCls++;}
 332   }
 333   if (sumCls>0) {
 334     qfactor = sumQ*1.0/sumCls;
 335     shfrac = sumSha*1.0/sumCls;
 336   }
 337   return qfactor;
 338 }
 339
 340 // Rotate second track for mixed pairs
 341 double AliTwoTrackRes::RotTr2Phi() {
 342   double rot = TVector2::Phi_mpi_pi(fP2.Phi()+TMath::Pi());
 343   fTpcEnt2.SetPhi(TVector2::Phi_mpi_pi(fTpcEnt2.Phi()+TMath::Pi()));
 344   return rot;}
 345
 346 // Fill NTuples
 347 void AliTwoTrackRes::FillNTuple1(double minsep, double sep, double corr,
 348                                  double qf, int ns1, int ns2) {
 349   fNTuple1->Fill(fP1.Pt(),fP1.Eta(),fP1.Phi(),ns1,fP2.Pt(),fP2.Eta(),
 350                  fP2.Phi(),ns2,Qinv(),minsep,sep,corr,qf);}
 351 void AliTwoTrackRes::FillNTuple2(double minsep, double sep, double corr,
 352                                  double qf, int ns1, int ns2) {
 353   fNTuple2->Fill(fPb1.Pt(),fPb1.Eta(),fPb1.Phi(),ns1,fPb2.Pt(),fPb2.Eta(),
 354                  fPb2.Phi(),ns2,Qinv(),minsep,sep,corr,qf);}
 355
 356 //______________________________________________________________________________
 357 // EOF
 358
 359