PWG1/TPC/AliPerformanceMatch.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceMatch class. It checks the track matching
   3 // quality (residuals, pulls) bewteen tracking detectors TPC-ITS and TPC-TRD
   4 // with TPC as the reference detector. In addition, the ITS and TRD track
   5 // reconstruction efficiency with respect to TPC is calculated.
   6 //
   7 // The matchinig quality parameters are stored in the THnSparse histograms.
   8 // The analysis of these histograms to extract reference information is done in
   9 // the Analyse() function (post-processing).
  10 //
  11 // The histograms with reference information can be exported to the ROOT folder.
  12 //
  13 // Author: J.Otwinowski 17/10/2009
  14 //------------------------------------------------------------------------------
  15
  16 /*
  17
  18   // after running comparison task, read the file, and get component
  19   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  20   LoadMyLibs();
  21
  22   TFile f("Output.root");
  23   AliPerformanceMatch * compObj = (AliPerformanceMatch*)coutput->FindObject("AliPerformanceMatchTPCTRD");
  24
  25   // analyse comparison data
  26   compObj->Analyse();
  27
  28   // the output histograms/graphs will be stored in the folder "folderMatch"
  29   compObj->GetAnalysisFolder()->ls("*");
  30
  31   // user can save whole comparison object (or only folder with anlysed histograms)
  32   // in the seperate output file (e.g.)
  33   TFile fout("Analysed_Match.root","recreate");
  34   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  35   fout.Close();
  36
  37 */
  38
  39 #include "TCanvas.h"
  40 #include "TH1.h"
  41 #include "TH2.h"
  42 #include "TAxis.h"
  43 #include "TF1.h"
  44
  45 #include "AliPerformanceMatch.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDVertex.h"
  48 #include "AliESDtrack.h"
  49 #include "AliESDfriendTrack.h"
  50 #include "AliESDfriend.h"
  51 #include "AliLog.h"
  52 #include "AliMCEvent.h"
  53 #include "AliMCParticle.h"
  54 #include "AliHeader.h"
  55 #include "AliGenEventHeader.h"
  56 #include "AliStack.h"
  57 #include "AliMCInfoCuts.h"
  58 #include "AliRecInfoCuts.h"
  59 #include "AliTracker.h"
  60 #include "AliTRDtrackV1.h"
  61 #include "AliTreeDraw.h"
  62
  63 using namespace std;
  64
  65 ClassImp(AliPerformanceMatch)
  66
  67 //_____________________________________________________________________________
  68 AliPerformanceMatch::AliPerformanceMatch():
  69   AliPerformanceObject("AliPerformanceMatch"),
  70   fResolHisto(0),
  71   fPullHisto(0),
  72   fTrackingEffHisto(0),
  73
  74   // Cuts
  75   fCutsRC(0),
  76   fCutsMC(0),
  77
  78   // histogram folder
  79   fAnalysisFolder(0)
  80 {
  81   Init();
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliPerformanceMatch::AliPerformanceMatch(Char_t* name="AliPerformanceMatch", Char_t* title="AliPerformanceMatch",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  86   AliPerformanceObject(name,title),
  87   fResolHisto(0),
  88   fPullHisto(0),
  89   fTrackingEffHisto(0),
  90
  91   // Cuts
  92   fCutsRC(0),
  93   fCutsMC(0),
  94
  95   // histogram folder
  96   fAnalysisFolder(0)
  97 {
  98   // named constructor
  99   //
 100   SetAnalysisMode(analysisMode);
 101   SetHptGenerator(hptGenerator);
 102
 103   Init();
 104 }
 105
 106 //_____________________________________________________________________________
 107 AliPerformanceMatch::~AliPerformanceMatch()
 108 {
 109   // destructor
 110
 111   if(fResolHisto) delete fResolHisto; fResolHisto=0;
 112   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;
 113   if(fTrackingEffHisto) delete fTrackingEffHisto; fTrackingEffHisto = 0x0;
 114
 115   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 116 }
 117
 118 //_____________________________________________________________________________
 119 void AliPerformanceMatch::Init(){
 120
 121   //
 122   // Make performance histogrms
 123   //
 124
 125   // set pt bins
 126   Int_t nPtBins = 50;
 127   Double_t ptMin = 1.e-2, ptMax = 10.;
 128
 129   Double_t *binsPt = 0;
 130   if (IsHptGenerator())  {
 131     nPtBins = 100; ptMax = 100.;
 132     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 133   } else {
 134     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 135   }
 136
 137   Double_t yMin = -0.02, yMax = 0.02;
 138   Double_t zMin = -12.0, zMax = 12.0;
 139   if(GetAnalysisMode() == 0 || GetAnalysisMode() == 1) {
 140     yMin = -100.; yMax = 100.;
 141     zMin = -100.; zMax = 100.;
 142   }
 143
 144   // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt:isRec
 145   Int_t binsResolHisto[11]={100,100,100,100,100,25,50,90,30,nPtBins,2};
 146   Double_t minResolHisto[11]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin, 0., -1.5, ptMin,0};
 147   Double_t maxResolHisto[11]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax,2};
 148
 149   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt:isRec",11,binsResolHisto,minResolHisto,maxResolHisto);
 150   fResolHisto->SetBinEdges(9,binsPt);
 151
 152   fResolHisto->GetAxis(0)->SetTitle("y-y_{ref} (cm)");
 153   fResolHisto->GetAxis(1)->SetTitle("z-z_{ref} (cm)");
 154   fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{ref} (rad)");
 155   fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{ref} (rad)");
 156   fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tref}-1)");
 157   fResolHisto->GetAxis(5)->SetTitle("y_{ref} (cm)");
 158   fResolHisto->GetAxis(6)->SetTitle("z_{ref} (cm)");
 159   fResolHisto->GetAxis(7)->SetTitle("#phi_{ref} (rad)");
 160   fResolHisto->GetAxis(8)->SetTitle("#eta_{ref}");
 161   fResolHisto->GetAxis(9)->SetTitle("p_{Tref} (GeV/c)");
 162   fResolHisto->GetAxis(10)->SetTitle("isReconstructed");
 163   fResolHisto->Sumw2();
 164
 165   //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt:isRec
 166   Int_t binsPullHisto[11]={100,100,100,100,100,50,50,50,50,nPtBins,2};
 167   Double_t minPullHisto[11]={-5.,-5.,-5.,-5.,-5., yMin, zMin,-1.,-2.0, ptMin,0};
 168   Double_t maxPullHisto[11]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax,2};
 169   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt:isRec",11,binsPullHisto,minPullHisto,maxPullHisto);
 170
 171   fPullHisto->GetAxis(0)->SetTitle("(y-y_{ref})/#sigma");
 172   fPullHisto->GetAxis(1)->SetTitle("(z-z_{ref})/#sigma");
 173   fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{ref})/#sigma");
 174   fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{ref})/#sigma");
 175   fPullHisto->GetAxis(4)->SetTitle("(p_{Tref}/p_{T}-1)/#sigma");
 176   fPullHisto->GetAxis(5)->SetTitle("y_{ref} (cm)");
 177   fPullHisto->GetAxis(6)->SetTitle("z_{ref} (cm)");
 178   fPullHisto->GetAxis(7)->SetTitle("sin#phi_{ref}");
 179   fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{ref}");
 180   fPullHisto->GetAxis(9)->SetTitle("1/p_{Tref} (GeV/c)^{-1}");
 181   fPullHisto->GetAxis(10)->SetTitle("isReconstructed");
 182   fPullHisto->Sumw2();
 183
 184   // -> has match:y:z:snp:tgl:phi:pt:ITSclusters
 185   Int_t binsTrackingEffHisto[8]    = { 2,    50, 100, 50, 50, 90,           100,  7   };
 186   Double_t minTrackingEffHisto[8]  = {-0.5, -25, -50, -1, -2, 0.,            0,   -0.5 };
 187   Double_t maxTrackingEffHisto[8]  = { 1.5,  25,  50,  1,  2, 2*TMath::Pi(), 20,   6.5 };
 188
 189   fTrackingEffHisto = new THnSparseF("fTrackingEffHisto","has match:y:z:snp:tgl:phi:pt:ITSclusters",8,binsTrackingEffHisto,minTrackingEffHisto,maxTrackingEffHisto);
 190   fTrackingEffHisto->GetAxis(0)->SetTitle("IsMatching");
 191   fTrackingEffHisto->GetAxis(1)->SetTitle("local y (cm)");
 192   fTrackingEffHisto->GetAxis(2)->SetTitle("z (cm)");
 193   fTrackingEffHisto->GetAxis(3)->SetTitle("sin(#phi)");
 194   fTrackingEffHisto->GetAxis(4)->SetTitle("tan(#lambda)");
 195   fTrackingEffHisto->GetAxis(5)->SetTitle("phi (rad)");
 196   fTrackingEffHisto->GetAxis(6)->SetTitle("p_{T}");
 197   fTrackingEffHisto->GetAxis(7)->SetTitle("number of ITS clusters");
 198   fTrackingEffHisto->Sumw2();
 199
 200   // Init cuts
 201   if(!fCutsMC)
 202     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 203   if(!fCutsRC)
 204     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 205
 206   // init folder
 207   fAnalysisFolder = CreateFolder("folderMatch","Analysis Matching Folder");
 208 }
 209
 210 //_____________________________________________________________________________
 211 void AliPerformanceMatch::ProcessITSTPC(Int_t iTrack, AliESDEvent *const esdEvent, AliStack* /*const stack*/, AliESDtrack *const esdTrack, AliESDfriendTrack *const esdFriendTrack)
 212 {
 213   //
 214   // addition to standard analysis - check if ITS stand-alone tracks have a match in the TPC
 215   // Origin: A. Kalwait
 216   // Modified: J. Otwinowski
 217   if(!esdEvent) return;
 218   if(!esdTrack) return;
 219   if(!esdFriendTrack) return;
 220
 221   // ITS stand alone tracks with SPD layers
 222   if(!(esdTrack->GetStatus() & AliESDtrack::kITSpureSA)) return;
 223   if(!(esdTrack->GetStatus() & AliESDtrack::kITSrefit)) return;
 224   if(esdTrack->GetNcls(0)<fCutsRC->GetMinNClustersITS()) return;
 225   if(!esdTrack->HasPointOnITSLayer(0) && !esdTrack->HasPointOnITSLayer(1)) return;
 226
 227   Bool_t hasMatch = kFALSE;
 228     for (Int_t jTrack = 0; jTrack < esdEvent->GetNumberOfTracks(); jTrack++) {
 229       // loop for all those tracks and check if the corresponding TPC track is found
 230       if (jTrack==iTrack) continue;
 231       AliESDtrack *trackTPC = esdEvent->GetTrack(jTrack);
 232       if (!trackTPC) continue;
 233       if (!trackTPC->GetTPCInnerParam()) continue;
 234
 235       // TPC nClust/track after first tracking pass
 236       if(trackTPC->GetTPCNclsIter1()<fCutsRC->GetMinNClustersTPC()) continue;
 237
 238       AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(trackTPC->GetTPCInnerParam()));
 239       if(!innerTPC) continue;
 240
 241       Double_t x[3]; trackTPC->GetXYZ(x);
 242       Double_t b[3]; AliTracker::GetBxByBz(x,b);
 243       Double_t dca[2],cov[3];
 244       Bool_t isTPCOK = innerTPC->PropagateToDCABxByBz(esdEvent->GetPrimaryVertexSPD(),b,kVeryBig,dca,cov);
 245       if(!isTPCOK) {
 246         if(innerTPC) delete innerTPC;
 247         continue;
 248       }
 249
 250       //
 251       // select primaries
 252       //
 253       Double_t dcaToVertex = -1;
 254       if( fCutsRC->GetDCAToVertex2D() )
 255       {
 256         dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()+dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 257       }
 258       if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) {
 259         delete innerTPC; continue; }
 260       if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) {
 261         delete innerTPC; continue; }
 262       if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) {
 263         delete innerTPC; continue; }
 264
 265       // rough checking if they match
 266       /*
 267       printf("+++++++++++++++++++++++++++++++++++++++++++++++ ");
 268       printf("esdTrack->GetY() %e, esdTrack->GetSnp() %e, esdTrack->GetTgl() %e \n",
 269               esdTrack->GetY(), esdTrack->GetSnp(), esdTrack->GetTgl());
 270       printf("innerTPC->GetY() %e, innerTPC->GetSnp() %e, innerTPC->GetTgl() %e \n",
 271               innerTPC->GetY() , innerTPC->GetSnp() , innerTPC->GetTgl());
 272       */
 273       if (TMath::Abs(esdTrack->GetY() - innerTPC->GetY()) > 3) {
 274         delete innerTPC; continue;
 275       }
 276       if (TMath::Abs(esdTrack->GetSnp() - innerTPC->GetSnp()) > 0.2) {
 277         delete innerTPC; continue;
 278       }
 279       if (TMath::Abs(esdTrack->GetTgl() - innerTPC->GetTgl()) > 0.2) {
 280         delete innerTPC; continue;
 281       }
 282
 283       hasMatch = kTRUE;
 284       if(innerTPC) delete innerTPC;
 285     }
 286     //has match:y:z:snp:tgl:phi:pt:ITSclusters
 287     Double_t vecTrackingEff[8] = { hasMatch,esdTrack->GetY(),esdTrack->GetZ(),esdTrack->GetSnp(),esdTrack->GetTgl(),esdTrack->Phi(), esdTrack->Pt(),esdTrack->GetITSclusters(0) };
 288     fTrackingEffHisto->Fill(vecTrackingEff);
 289
 290 }
 291
 292 //_____________________________________________________________________________
 293 void AliPerformanceMatch::ProcessTPCITS(AliStack* /*const stack*/, AliESDtrack *const esdTrack, AliESDfriendTrack *const esdFriendTrack)
 294 {
 295   //
 296   // Match TPC and ITS min-bias tracks
 297   // at radius between detectors
 298   //
 299   if(!esdTrack) return;
 300   if(!esdFriendTrack) return;
 301
 302   //
 303   // Propagate tracks to the radius between TPC-ITS
 304   // using B-field and material budget
 305   //
 306   Double_t radius = fCutsRC->GetTPCITSMatchingRadius();
 307   Double_t mass = esdTrack->GetMass();
 308   Double_t step=1.0; // cm
 309
 310   //
 311   // Propagate TPCinner (reference detector)
 312   //
 313   Bool_t isTPCOK=kFALSE;
 314   AliExternalTrackParam *innerTPC=NULL;
 315
 316   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 317   esdTrack->GetImpactParametersTPC(dca,cov);
 318
 319   //
 320   // select primaries
 321   //
 322   Double_t dcaToVertex = -1;
 323   if( fCutsRC->GetDCAToVertex2D() )
 324   {
 325       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 326   }
 327   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 328   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 329   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 330
 331   if( (esdTrack->GetTPCNclsIter1()>fCutsRC->GetMinNClustersTPC()) &&
 332       (esdTrack->GetTPCInnerParam()) &&
 333       (innerTPC=new AliExternalTrackParam(*(esdTrack->GetTPCInnerParam()))))
 334   {
 335      isTPCOK = AliTracker::PropagateTrackToBxByBz(innerTPC,radius,mass,step,kTRUE);
 336   }
 337   if(!isTPCOK) {
 338    if(innerTPC) delete innerTPC;
 339    return;
 340   }
 341
 342   //
 343   // Propagate ITSouter
 344   //
 345   Bool_t isITSOK=kFALSE;
 346   AliExternalTrackParam *outerITS=NULL;
 347
 348   if( (esdTrack->GetNcls(0)>fCutsRC->GetMinNClustersITS()) &&
 349       (esdFriendTrack->GetITSOut()) &&
 350       (outerITS=new AliExternalTrackParam(*(esdFriendTrack->GetITSOut()))))
 351   {
 352     isITSOK = AliTracker::PropagateTrackToBxByBz(outerITS,radius,mass,step,kTRUE);
 353   }
 354
 355   //
 356   // Fill histograms (TPC reference detector)
 357   //
 358   if(isTPCOK)
 359     FillHistograms(innerTPC,outerITS,isITSOK);
 360
 361   if(outerITS) delete outerITS;
 362   if(innerTPC) delete innerTPC;
 363 }
 364
 365 //_____________________________________________________________________________
 366 void AliPerformanceMatch::ProcessTPCTRD(AliStack* /*const stack*/, AliESDtrack *const esdTrack, AliESDfriendTrack *const esdFriendTrack)
 367 {
 368   //
 369   // Match TPC and TRD min-bias tracks
 370   // at radius between detectors. TPC is the reference detector.
 371   //
 372   if(!esdTrack) return;
 373   if(!esdFriendTrack) return;
 374
 375   //
 376   // Propagate tracks to the radius between TPC-TRD
 377   // using B-field and material budget
 378   //
 379   Double_t radius = fCutsRC->GetTPCTRDMatchingRadius();
 380   Double_t mass = esdTrack->GetMass();
 381   Double_t step=1.0; // cm
 382
 383   //
 384   // Propagate TPCouter (reference detector)
 385   //
 386   Bool_t isTPCOK=kFALSE;
 387   AliExternalTrackParam *outerTPC=NULL;
 388
 389   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 390   esdTrack->GetImpactParametersTPC(dca,cov);
 391
 392   //
 393   // select primaries
 394   //
 395   Double_t dcaToVertex = -1;
 396   if( fCutsRC->GetDCAToVertex2D() )
 397   {
 398       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 399   }
 400   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 401   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 402   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 403
 404   if( (esdTrack->GetTPCNclsIter1()>fCutsRC->GetMinNClustersTPC()) &&
 405       (esdFriendTrack->GetTPCOut()) &&
 406       (outerTPC=new AliExternalTrackParam(*(esdFriendTrack->GetTPCOut()))))
 407   {
 408      isTPCOK = AliTracker::PropagateTrackToBxByBz(outerTPC,radius,mass,step,kTRUE);
 409   }
 410   if(!isTPCOK)  {
 411     if(outerTPC) delete outerTPC;
 412     return;
 413   }
 414
 415   //
 416   // Propagate TRDinner
 417   //
 418   Bool_t isTRDOK = kFALSE;
 419   AliExternalTrackParam *innerTRD=NULL;
 420
 421   // get TRD track
 422   AliTRDtrackV1 *trdTrack=NULL; //esdFriendTrack = fESDfriend->GetTrack(itrk);
 423   TObject *calObject=NULL;
 424   Int_t icalib = 0;
 425   while((calObject = esdFriendTrack->GetCalibObject(icalib++))) {
 426     if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
 427     if(!(trdTrack = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
 428   }
 429
 430   if( (trdTrack) &&
 431       (trdTrack->GetNumberOfTracklets()>fCutsRC->GetMinNTrackletsTRD()) &&
 432       (trdTrack->GetTracklet(0)) &&
 433       (esdFriendTrack->GetTRDIn()) &&
 434       (innerTRD = new AliExternalTrackParam(*(esdFriendTrack->GetTRDIn()))))
 435   {
 436     isTRDOK = AliTracker::PropagateTrackToBxByBz(innerTRD,radius,mass,step,kTRUE);
 437   }
 438
 439   //
 440   // Fill histograms (TPC reference detector)
 441   //
 442   if(isTPCOK)
 443     FillHistograms(outerTPC,innerTRD,isTRDOK);
 444
 445   if(outerTPC) delete outerTPC;
 446   if(innerTRD) delete innerTRD;
 447 }
 448
 449 //_____________________________________________________________________________
 450 void AliPerformanceMatch::FillHistograms(AliExternalTrackParam *const refParam, AliExternalTrackParam *const param, Bool_t isRec)
 451 {
 452   //
 453   // fill performance histograms
 454   // (TPC always as reference)
 455   //
 456   if(!refParam) return;
 457
 458   //
 459   // Deltas (dy,dz,dphi,dtheta,dpt)
 460   //
 461   Float_t delta[5] = {0};
 462   if(param && isRec) {
 463     delta[0] = param->GetY()-refParam->GetY();
 464     delta[1] = param->GetZ()-refParam->GetZ();
 465     delta[2] = TMath::ATan2(param->Py(),param->Px())-TMath::ATan2(refParam->Py(),refParam->Px());
 466     delta[3] = TMath::ATan2(param->Pz(),param->Pt())-TMath::ATan2(refParam->Pz(),refParam->Pt());
 467     if(refParam->Pt()) delta[4] = (param->Pt()-refParam->Pt())/refParam->Pt();
 468   }
 469   //
 470   // Pulls (y,z,snp,tanl,1/pt)
 471   //
 472   Float_t sigma[5] = {0};
 473   Float_t pull[5] = {0};
 474   if(param && isRec) {
 475     sigma[0] = TMath::Sqrt(param->GetSigmaY2()+refParam->GetSigmaY2());
 476     sigma[1] = TMath::Sqrt(param->GetSigmaZ2()+refParam->GetSigmaZ2());
 477     sigma[2] = TMath::Sqrt(param->GetSigmaSnp2()+refParam->GetSigmaSnp2());
 478     sigma[3] = TMath::Sqrt(param->GetSigmaTgl2()+refParam->GetSigmaTgl2());
 479     sigma[4] = TMath::Sqrt(param->GetSigma1Pt2()+refParam->GetSigma1Pt2());
 480     if(sigma[0]) pull[0] = delta[0] / sigma[0];
 481     if(sigma[1]) pull[1] = delta[1] / sigma[1];
 482     if(sigma[2]) pull[2] = (param->GetSnp()-refParam->GetSnp()) / sigma[2];
 483     if(sigma[3]) pull[3] = (param->GetTgl()-refParam->GetTgl()) / sigma[3];
 484     if(sigma[4]) pull[4] = (param->OneOverPt()-refParam->OneOverPt()) / sigma[4];
 485   }
 486
 487   // Fill histograms
 488   Double_t vResolHisto[11] = {delta[0],delta[1],delta[2],delta[3],delta[4],refParam->GetY(),refParam->GetZ(),refParam->Phi(),refParam->Eta(),refParam->Pt(),isRec};
 489   fResolHisto->Fill(vResolHisto);
 490
 491   Double_t vPullHisto[11] = {pull[0],pull[1],pull[2],pull[3],pull[4],refParam->GetY(),refParam->GetZ(),refParam->GetSnp(),refParam->GetTgl(),refParam->OneOverPt(),isRec};
 492   fPullHisto->Fill(vPullHisto);
 493 }
 494
 495 //_____________________________________________________________________________
 496 void AliPerformanceMatch::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 497 {
 498   // Process comparison information
 499   //
 500   if(!esdEvent)
 501   {
 502     Error("Exec","esdEvent not available");
 503     return;
 504   }
 505   AliHeader* header = 0;
 506   AliGenEventHeader* genHeader = 0;
 507   AliStack* stack = 0;
 508   TArrayF vtxMC(3);
 509
 510   if(bUseMC)
 511   {
 512     if(!mcEvent) {
 513       Error("Exec","mcEvent not available");
 514       return;
 515     }
 516     // get MC event header
 517     header = mcEvent->Header();
 518     if (!header) {
 519       Error("Exec","Header not available");
 520       return;
 521     }
 522     // MC particle stack
 523     stack = mcEvent->Stack();
 524     if (!stack) {
 525       Error("Exec","Stack not available");
 526       return;
 527     }
 528     // get MC vertex
 529     genHeader = header->GenEventHeader();
 530     if (!genHeader) {
 531       Error("Exec","Could not retrieve genHeader from Header");
 532       return;
 533     }
 534     genHeader->PrimaryVertex(vtxMC);
 535   }
 536
 537   // use ESD friends
 538   if(bUseESDfriend) {
 539     if(!esdFriend) {
 540       Error("Exec","esdFriend not available");
 541       return;
 542     }
 543   }
 544
 545   // trigger
 546   if(!bUseMC && GetTriggerClass()) {
 547     Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
 548     if(!isEventTriggered) return;
 549   }
 550
 551   // get TPC event vertex
 552   const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
 553   if(vtxESD && (vtxESD->GetStatus()<=0)) return;
 554
 555   //  Process events
 556   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 557   {
 558     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 559     if(!track) continue;
 560
 561     AliESDfriendTrack *friendTrack=0;
 562     if(bUseESDfriend) {
 563       friendTrack=esdFriend->GetTrack(iTrack);
 564       if(!friendTrack) continue;
 565     }
 566
 567     if(GetAnalysisMode() == 0) ProcessTPCITS(stack,track,friendTrack);
 568     else if(GetAnalysisMode() == 1) ProcessTPCTRD(stack,track,friendTrack);
 569     else if(GetAnalysisMode() == 2) ProcessITSTPC(iTrack,esdEvent,stack,track,friendTrack);
 570     else {
 571       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 572       return;
 573     }
 574   }
 575 }
 576
 577 //_____________________________________________________________________________
 578 TH1F* AliPerformanceMatch::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
 579   // Create resolution histograms
 580
 581   TH1F *hisr, *hism;
 582   if (!gPad) new TCanvas;
 583   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
 584   if (type) return hism;
 585   else
 586     return hisr;
 587 }
 588
 589 //_____________________________________________________________________________
 590 void AliPerformanceMatch::Analyse() {
 591   // Analyse comparison information and store output histograms
 592   // in the folder "folderMatch"
 593   //
 594   TH1::AddDirectory(kFALSE);
 595   TH1F *h=0;
 596   TH1F *h2=0;
 597   TH2F *h2D=0;
 598   TObjArray *aFolderObj = new TObjArray;
 599
 600   // write results in the folder
 601   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 602   c->cd();
 603
 604   char name[256];
 605   char title[256];
 606
 607   if(GetAnalysisMode()==0 || GetAnalysisMode()==1) {
 608
 609   fResolHisto->GetAxis(10)->SetRangeUser(1.0,2.0); // only reconstructed
 610   fPullHisto->GetAxis(10)->SetRangeUser(1.0,2.0);  // only reconstructed
 611   for(Int_t i=0; i<5; i++)
 612   {
 613     for(Int_t j=5; j<10; j++)
 614     {
 615       //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
 616       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
 617       else fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);
 618       fResolHisto->GetAxis(9)->SetRangeUser(0.1,100.); // pt threshold
 619
 620       h2D = (TH2F*)fResolHisto->Projection(i,j);
 621
 622       h = AliPerformanceMatch::MakeResol(h2D,1,0,100);
 623       sprintf(name,"h_res_%d_vs_%d",i,j);
 624       h->SetName(name);
 625
 626       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 627       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
 628       h->GetYaxis()->SetTitle(title);
 629       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 630       h->SetTitle(title);
 631
 632       //if(j==9) h->SetBit(TH1::kLogX);
 633       aFolderObj->Add(h);
 634
 635       h = AliPerformanceMatch::MakeResol(h2D,1,1,100);
 636       //h = (TH1F*)arr->At(1);
 637       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
 638       h->SetName(name);
 639
 640       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 641       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
 642       h->GetYaxis()->SetTitle(title);
 643
 644       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 645       h->SetTitle(title);
 646
 647       if(j==9) h->SetBit(TH1::kLogX);
 648       aFolderObj->Add(h);
 649
 650       //
 651       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
 652       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.49); // eta window
 653       fPullHisto->GetAxis(9)->SetRangeUser(0.1,100.);  // pt threshold
 654
 655       h2D = (TH2F*)fPullHisto->Projection(i,j);
 656
 657       h = AliPerformanceMatch::MakeResol(h2D,1,0,100);
 658       sprintf(name,"h_pull_%d_vs_%d",i,j);
 659       h->SetName(name);
 660
 661       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 662       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
 663       h->GetYaxis()->SetTitle(title);
 664       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 665       h->SetTitle(title);
 666
 667       if(j==9) h->SetBit(TH1::kLogX);
 668       aFolderObj->Add(h);
 669
 670       h = AliPerformanceMatch::MakeResol(h2D,1,1,100);
 671       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
 672       h->SetName(name);
 673
 674       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 675       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
 676       h->GetYaxis()->SetTitle(title);
 677       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 678       h->SetTitle(title);
 679
 680       //if(j==9) h->SetBit(TH1::kLogX);
 681       aFolderObj->Add(h);
 682     }
 683   }
 684
 685   //
 686   // Efficiency plots
 687   //
 688   for(Int_t i=5; i<10; i++)
 689   {
 690     if(i!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
 691     else fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);
 692     fResolHisto->GetAxis(9)->SetRangeUser(0.1,100.); // pt threshold
 693
 694     fResolHisto->GetAxis(10)->SetRange(1,fResolHisto->GetAxis(10)->GetNbins()); // all
 695     h = (TH1F*)fResolHisto->Projection(i);
 696
 697     fResolHisto->GetAxis(10)->SetRange(2,2); // only reconstructed
 698     h2 = (TH1F*)fResolHisto->Projection(i);
 699
 700     TH1F* h2c = (TH1F*)h2->Clone();
 701     h2c->Divide(h2,h,1,1,"B");
 702
 703     sprintf(name,"h_eff_%d",i);
 704     h2c->SetName(name);
 705
 706     h2c->GetXaxis()->SetTitle(h2c->GetXaxis()->GetTitle());
 707     h2c->GetYaxis()->SetTitle("efficiency");
 708     h2c->SetTitle("matching effciency");
 709
 710     aFolderObj->Add(h2c);
 711   }
 712
 713   }
 714
 715   //
 716   // TPC efficiency wrt ITS
 717   //
 718   if(GetAnalysisMode()==2) {
 719
 720     h = (TH1F*)fTrackingEffHisto->Projection(0);
 721     aFolderObj->Add(h);
 722
 723     for(Int_t i=1; i<7; i++)
 724     {
 725       //
 726       //
 727       // calculate efficiency
 728       //
 729
 730       // all ITS standalone tracks
 731       fTrackingEffHisto->GetAxis(0)->SetRange(1,fTrackingEffHisto->GetAxis(0)->GetNbins());
 732       h = (TH1F*)fTrackingEffHisto->Projection(i);
 733
 734       // TPC tracks which has matching with TPC
 735       fTrackingEffHisto->GetAxis(0)->SetRange(2,2);
 736       h2 = (TH1F*)fTrackingEffHisto->Projection(i);
 737
 738       TH1F* h2c = (TH1F*)h2->Clone();
 739       h2c->Divide(h2,h,1,1,"B");
 740
 741       sprintf(name,"h_TPC_eff_%d",i);
 742       h2c->SetName(name);
 743
 744       h2c->GetXaxis()->SetTitle(h2c->GetXaxis()->GetTitle());
 745       h2c->GetYaxis()->SetTitle("efficiency");
 746       h2c->SetTitle("TPC effciency wrt ITS");
 747
 748       aFolderObj->Add(h2c);
 749     }
 750
 751   }
 752
 753   // export objects to analysis folder
 754   fAnalysisFolder = ExportToFolder(aFolderObj);
 755
 756   // delete only TObjArray
 757   if(aFolderObj) delete aFolderObj;
 758 }
 759
 760 //_____________________________________________________________________________
 761 TFolder* AliPerformanceMatch::ExportToFolder(TObjArray * array)
 762 {
 763   // recreate folder avery time and export objects to new one
 764   //
 765   AliPerformanceMatch * comp=this;
 766   TFolder *folder = comp->GetAnalysisFolder();
 767
 768   TString name, title;
 769   TFolder *newFolder = 0;
 770   Int_t i = 0;
 771   Int_t size = array->GetSize();
 772
 773   if(folder) {
 774      // get name and title from old folder
 775      name = folder->GetName();
 776      title = folder->GetTitle();
 777
 778          // delete old one
 779      delete folder;
 780
 781          // create new one
 782      newFolder = CreateFolder(name.Data(),title.Data());
 783      newFolder->SetOwner();
 784
 785          // add objects to folder
 786      while(i < size) {
 787            newFolder->Add(array->At(i));
 788            i++;
 789          }
 790   }
 791
 792 return newFolder;
 793 }
 794
 795 //_____________________________________________________________________________
 796 Long64_t AliPerformanceMatch::Merge(TCollection* const list)
 797 {
 798   // Merge list of objects (needed by PROOF)
 799
 800   if (!list)
 801   return 0;
 802
 803   if (list->IsEmpty())
 804   return 1;
 805
 806   TIterator* iter = list->MakeIterator();
 807   TObject* obj = 0;
 808
 809   // collection of generated histograms
 810   Int_t count=0;
 811   while((obj = iter->Next()) != 0)
 812   {
 813   AliPerformanceMatch* entry = dynamic_cast<AliPerformanceMatch*>(obj);
 814   if (entry == 0) continue;
 815
 816   fResolHisto->Add(entry->fResolHisto);
 817   fPullHisto->Add(entry->fPullHisto);
 818   fTrackingEffHisto->Add(entry->fTrackingEffHisto);
 819
 820   count++;
 821   }
 822
 823 return count;
 824 }
 825
 826 //_____________________________________________________________________________
 827 TFolder* AliPerformanceMatch::CreateFolder(TString name,TString title) {
 828 // create folder for analysed histograms
 829 //
 830 TFolder *folder = 0;
 831   folder = new TFolder(name.Data(),title.Data());
 832
 833   return folder;
 834 }