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