PWG1/TPC/AliPerformanceTPC.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceTPC class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms/graphs) are stored in the folder which is
   8 // a data member of AliPerformanceTPC.
   9 //
  10 // Author: J.Otwinowski 04/02/2008
  11 // Changes by M.Knichel 15/10/2010
  12 //------------------------------------------------------------------------------
  13
  14 /*
  15
  16   // after running comparison task, read the file, and get component
  17   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  18   LoadMyLibs();
  19
  20   TFile f("Output.root");
  21   AliPerformanceTPC * compObj = (AliPerformanceTPC*)coutput->FindObject("AliPerformanceTPC");
  22
  23   // analyse comparison data
  24   compObj->Analyse();
  25
  26   // the output histograms/graphs will be stored in the folder "folderTPC"
  27   compObj->GetAnalysisFolder()->ls("*");
  28
  29   // user can save whole comparison object (or only folder with anlysed histograms)
  30   // in the seperate output file (e.g.)
  31   TFile fout("Analysed_TPC.root","recreate");
  32   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  33   fout.Close();
  34
  35 */
  36
  37 #include "TCanvas.h"
  38 #include "TH1.h"
  39 #include "TH2.h"
  40 #include "TH3.h"
  41 #include "TAxis.h"
  42 #include "TPostScript.h"
  43 #include "TString.h"
  44 #include "TUUID.h"
  45 #include "TTree.h"
  46 #include "TChain.h"
  47 #include "AliTPCPerformanceSummary.h"
  48 #include "TSystem.h"
  49
  50 #include "AliPerformanceTPC.h"
  51 #include "AliESDEvent.h"
  52 #include "AliESDVertex.h"
  53 #include "AliESDtrack.h"
  54 #include "AliLog.h"
  55 #include "AliMCEvent.h"
  56 #include "AliHeader.h"
  57 #include "AliGenEventHeader.h"
  58 #include "AliStack.h"
  59 #include "AliMCInfoCuts.h"
  60 #include "AliRecInfoCuts.h"
  61 #include "AliTracker.h"
  62 #include "AliTreeDraw.h"
  63 #include "AliTPCTransform.h"
  64 #include "AliTPCseed.h"
  65 #include "AliTPCcalibDB.h"
  66 #include "AliESDfriend.h"
  67 #include "AliESDfriendTrack.h"
  68 #include "AliTPCclusterMI.h"
  69
  70 using namespace std;
  71
  72 ClassImp(AliPerformanceTPC)
  73
  74 Bool_t AliPerformanceTPC::fgMergeTHnSparse = kFALSE;
  75
  76
  77 //_____________________________________________________________________________
  78 /*
  79 AliPerformanceTPC::AliPerformanceTPC():
  80   AliPerformanceObject("AliPerformanceTPC"),
  81   fTPCClustHisto(0),
  82   fTPCEventHisto(0),
  83   fTPCTrackHisto(0),
  84   fFolderObj(0),
  85
  86   // Cuts
  87   fCutsRC(0),
  88   fCutsMC(0),
  89
  90   // histogram folder
  91   fAnalysisFolder(0),
  92
  93   fUseHLT(kFALSE)
  94
  95 {
  96   Init();
  97 }
  98 */
  99
 100 //_____________________________________________________________________________
 101 AliPerformanceTPC::AliPerformanceTPC(Char_t* name, Char_t* title,Int_t analysisMode,Bool_t hptGenerator, Int_t run, Bool_t highMult):
 102   AliPerformanceObject(name,title,run,highMult),
 103   fTPCClustHisto(0),
 104   fTPCEventHisto(0),
 105   fTPCTrackHisto(0),
 106   fFolderObj(0),
 107
 108   // Cuts
 109   fCutsRC(0),
 110   fCutsMC(0),
 111
 112   // histogram folder
 113   fAnalysisFolder(0),
 114
 115   fUseHLT(kFALSE)
 116
 117 {
 118   // named constructor
 119   //
 120   SetAnalysisMode(analysisMode);
 121   SetHptGenerator(hptGenerator);
 122
 123   Init();
 124 }
 125
 126
 127 //_____________________________________________________________________________
 128 AliPerformanceTPC::~AliPerformanceTPC()
 129 {
 130   // destructor
 131
 132   if(fTPCClustHisto) delete fTPCClustHisto; fTPCClustHisto=0;
 133   if(fTPCEventHisto) delete fTPCEventHisto; fTPCEventHisto=0;
 134   if(fTPCTrackHisto) delete fTPCTrackHisto; fTPCTrackHisto=0;
 135   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 136   if(fFolderObj) delete fFolderObj; fFolderObj=0;
 137 }
 138
 139
 140 //_____________________________________________________________________________
 141 void AliPerformanceTPC::Init()
 142 {
 143   //
 144   // histogram bining
 145   //
 146
 147   // set pt bins
 148   Int_t nPtBins = 50;
 149   Double_t ptMin = 1.e-2, ptMax = 20.;
 150
 151   Double_t *binsPt = 0;
 152   if (IsHptGenerator())  {
 153     nPtBins = 100; ptMax = 100.;
 154     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 155   } else {
 156     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 157   }
 158
 159   const Int_t  nCOverPtBins = 80;
 160   Double_t coverptMin = -10, coverptMax = 10;
 161   Double_t *binsCOverPtP = 0;
 162   Double_t *binsCOverPt = new Double_t[nCOverPtBins+1];
 163   binsCOverPtP = CreateLogAxis(nCOverPtBins/2,0.04,coverptMax-0.04);
 164   for(Int_t i=0; i < nCOverPtBins/2; i++){
 165     binsCOverPt[nCOverPtBins - i] = binsCOverPtP[nCOverPtBins/2  - i];
 166     binsCOverPt[i] = 0 - binsCOverPtP[nCOverPtBins/2  - i];
 167  }
 168
 169   /*
 170   Int_t nPtBins = 31;
 171   Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
 172   Double_t ptMin = 0., ptMax = 10.;
 173
 174   if(IsHptGenerator() == kTRUE) {
 175     nPtBins = 100;
 176     ptMin = 0.; ptMax = 100.;
 177   }
 178   */
 179   //
 180
 181   //
 182   //padRow:phi:TPCSide
 183   Int_t binsTPCClustHisto[3] =   {160,  144,  2};
 184   Double_t minTPCClustHisto[3] = {0.,   0.,   0.};
 185   Double_t maxTPCClustHisto[3] = {160., 2.*TMath::Pi(), 2.};
 186
 187   fTPCClustHisto = new THnSparseF("fTPCClustHisto","padRow:phi:TPCSide",3,binsTPCClustHisto,minTPCClustHisto,maxTPCClustHisto);
 188   fTPCClustHisto->GetAxis(0)->SetTitle("padRow");
 189   fTPCClustHisto->GetAxis(1)->SetTitle("phi (rad)");
 190   fTPCClustHisto->GetAxis(2)->SetTitle("TPCSide");
 191   //fTPCClustHisto->Sumw2();
 192
 193   //padRow:phi:TPCSide:pad:detector:glZ
 194   /*
 195   Int_t binsTPCClustHisto[6] =   {160,  144,  2, 128, 72, 50};
 196   Double_t minTPCClustHisto[6] = {0.,   0.,   0., 0, 0, -250};
 197   Double_t maxTPCClustHisto[6] = {160., 2.*TMath::Pi(), 2., 128, 72,250};
 198
 199   fTPCClustHisto = new THnSparseF("fTPCClustHisto","padRow:phi:TPCSide:pad:detector:gZ",6,binsTPCClustHisto,minTPCClustHisto,maxTPCClustHisto);
 200   fTPCClustHisto->GetAxis(0)->SetTitle("padRow");
 201   fTPCClustHisto->GetAxis(1)->SetTitle("phi (rad)");
 202   fTPCClustHisto->GetAxis(2)->SetTitle("TPCSide");
 203   fTPCClustHisto->GetAxis(3)->SetTitle("pad");
 204   fTPCClustHisto->GetAxis(4)->SetTitle("detector");
 205   fTPCClustHisto->GetAxis(5)->SetTitle("glZ (cm)");
 206   //fTPCClustHisto->Sumw2();
 207   */
 208
 209   Float_t scaleVxy = 1.0;
 210   if(fAnalysisMode !=0) scaleVxy = 0.1;
 211
 212   Int_t maxMult;
 213   if (fHighMultiplicity) { maxMult = 4001; scaleVxy = 0.1;} else { maxMult = 151; }
 214   // Xv:Yv:Zv:mult:multP:multN:vertStatus
 215   Int_t binsTPCEventHisto[7]=  {100,  100,   100,  maxMult,  maxMult,  maxMult, 2   };
 216   Double_t minTPCEventHisto[7]={-10.*scaleVxy, -10.*scaleVxy, -30., -0.5,  -0.5,  -0.5, -0.5  };
 217   Double_t maxTPCEventHisto[7]={ 10.*scaleVxy,  10.*scaleVxy,  30.,  maxMult-0.5,  maxMult-0.5, maxMult-0.5, 1.5 };
 218
 219   fTPCEventHisto = new THnSparseF("fTPCEventHisto","Xv:Yv:Zv:mult:multP:multN:vertStatus",7,binsTPCEventHisto,minTPCEventHisto,maxTPCEventHisto);
 220   fTPCEventHisto->GetAxis(0)->SetTitle("Xv (cm)");
 221   fTPCEventHisto->GetAxis(1)->SetTitle("Yv (cm)");
 222   fTPCEventHisto->GetAxis(2)->SetTitle("Zv (cm)");
 223   fTPCEventHisto->GetAxis(3)->SetTitle("mult");
 224   fTPCEventHisto->GetAxis(4)->SetTitle("multP");
 225   fTPCEventHisto->GetAxis(5)->SetTitle("multN");
 226   fTPCEventHisto->GetAxis(6)->SetTitle("vertStatus");
 227   //fTPCEventHisto->Sumw2();
 228
 229   Float_t scaleDCA = 1.0;
 230   if(fAnalysisMode !=0) scaleDCA = 0.1;
 231   // nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:DCAr:DCAz:eta:phi:pt:charge:vertStatus
 232    //Int_t binsTPCTrackHisto[10]=  { 160,  20,  60,  30, 30,  30,   144,             nPtBins,   nCOverPtBins, 2 };
 233    //Double_t minTPCTrackHisto[10]={ 0.,   0.,  0., -3*scaleDCA, -3.*scaleDCA, -1.5, 0.,             ptMin,   coverptMin, -0.5 };
 234    //Double_t maxTPCTrackHisto[10]={ 160., 5., 1.2, 3*scaleDCA,  3.*scaleDCA,  1.5, 2.*TMath::Pi(), ptMax,    coverptMax,  1.5 };
 235    Int_t binsTPCTrackHisto[10]=  { 160,  20,  60,  30, 30,  30,   144,             nPtBins,   3, 2 };
 236    Double_t minTPCTrackHisto[10]={ 0.,   0.,  0., -3*scaleDCA, -3.*scaleDCA, -1.5, 0.,             ptMin,  -1.5, -0.5 };
 237    Double_t maxTPCTrackHisto[10]={ 160., 5., 1.2, 3*scaleDCA,  3.*scaleDCA,  1.5, 2.*TMath::Pi(), ptMax,    1.5,  1.5 };
 238
 239   fTPCTrackHisto = new THnSparseF("fTPCTrackHisto","nClust:chi2PerClust:nClust/nFindableClust:DCAr:DCAz:eta:phi:pt:charge/pt:vertStatus",10,binsTPCTrackHisto,minTPCTrackHisto,maxTPCTrackHisto);
 240   fTPCTrackHisto->SetBinEdges(7,binsPt);
 241   //fTPCTrackHisto->SetBinEdges(8,binsCOverPt);
 242
 243   fTPCTrackHisto->GetAxis(0)->SetTitle("nClust");
 244   fTPCTrackHisto->GetAxis(1)->SetTitle("chi2PerClust");
 245   fTPCTrackHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
 246   fTPCTrackHisto->GetAxis(3)->SetTitle("DCAr (cm)");
 247   fTPCTrackHisto->GetAxis(4)->SetTitle("DCAz (cm)");
 248   fTPCTrackHisto->GetAxis(5)->SetTitle("#eta");
 249   fTPCTrackHisto->GetAxis(6)->SetTitle("#phi (rad)");
 250   fTPCTrackHisto->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
 251   //fTPCTrackHisto->GetAxis(8)->SetTitle("charge/pt");
 252   fTPCTrackHisto->GetAxis(8)->SetTitle("charge");
 253   fTPCTrackHisto->GetAxis(9)->SetTitle("vertStatus");
 254   //fTPCTrackHisto->Sumw2();
 255
 256   // Init cuts
 257   if(!fCutsMC) {
 258     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 259   }
 260   if(!fCutsRC) {
 261     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 262   }
 263
 264   // init folder
 265   fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
 266
 267   delete []binsCOverPt;
 268 }
 269
 270
 271 //_____________________________________________________________________________
 272 void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent *const esdEvent, Bool_t vertStatus)
 273 {
 274 //
 275 // fill TPC QA info
 276 //
 277   if(!esdEvent) return;
 278   if(!esdTrack) return;
 279
 280   if( IsUseTrackVertex() )
 281   {
 282     // Relate TPC inner params to prim. vertex
 283     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 284     Double_t x[3]; esdTrack->GetXYZ(x);
 285     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 286     Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
 287     if(!isOK) return;
 288
 289     /*
 290       // JMT -- recaluclate DCA for HLT if not present
 291       if ( dca[0] == 0. && dca[1] == 0. ) {
 292         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 293       }
 294     */
 295   }
 296
 297   // Fill TPC only resolution comparison information
 298   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
 299   if(!track) return;
 300
 301   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 302   esdTrack->GetImpactParametersTPC(dca,cov);
 303
 304   Float_t q = esdTrack->Charge();
 305   Float_t pt = track->Pt();
 306   Float_t eta = track->Eta();
 307   Float_t phi = track->Phi();
 308   Int_t nClust = esdTrack->GetTPCclusters(0);
 309   Int_t nFindableClust = esdTrack->GetTPCNclsF();
 310
 311   Float_t chi2PerCluster = 0.;
 312   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
 313
 314   Float_t clustPerFindClust = 0.;
 315   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
 316
 317   Float_t qpt = 0;
 318   if( fabs(pt)>0 ) qpt = q/fabs(pt);
 319
 320   // filter out noise tracks
 321   if(esdTrack->GetTPCsignal() < 5) return;
 322
 323   //
 324   // select primaries
 325   //
 326   Double_t dcaToVertex = -1;
 327   if( fCutsRC->GetDCAToVertex2D() )
 328   {
 329       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 330   }
 331   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 332   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 333   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 334
 335   Double_t vTPCTrackHisto[10] = {nClust,chi2PerCluster,clustPerFindClust,dca[0],dca[1],eta,phi,pt,qpt,vertStatus};
 336   fTPCTrackHisto->Fill(vTPCTrackHisto);
 337
 338   //
 339   // Fill rec vs MC information
 340   //
 341   if(!stack) return;
 342
 343 }
 344
 345
 346 //_____________________________________________________________________________
 347 void AliPerformanceTPC::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent, Bool_t vertStatus)
 348 {
 349   // Fill comparison information (TPC+ITS)
 350   if(!esdTrack) return;
 351   if(!esdEvent) return;
 352
 353   if( IsUseTrackVertex() )
 354   {
 355     // Relate TPC inner params to prim. vertex
 356     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 357     Double_t x[3]; esdTrack->GetXYZ(x);
 358     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 359     Bool_t isOK = esdTrack->RelateToVertexBxByBz(vtxESD, b, kVeryBig);
 360     if(!isOK) return;
 361
 362     /*
 363       // JMT -- recaluclate DCA for HLT if not present
 364       if ( dca[0] == 0. && dca[1] == 0. ) {
 365         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 366       }
 367     */
 368   }
 369
 370   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 371   esdTrack->GetImpactParameters(dca,cov);
 372
 373   if ((esdTrack->GetStatus()&AliESDtrack::kITSrefit)==0) return; // ITS refit
 374   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 375   if ((esdTrack->HasPointOnITSLayer(0)==kFALSE)&&(esdTrack->HasPointOnITSLayer(1)==kFALSE)) return; // at least one SPD
 376   //if (esdTrack->GetITSclusters(0)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters
 377
 378   Float_t q = esdTrack->Charge();
 379   Float_t pt = esdTrack->Pt();
 380   Float_t eta = esdTrack->Eta();
 381   Float_t phi = esdTrack->Phi();
 382   Int_t nClust = esdTrack->GetTPCclusters(0);
 383   Int_t nFindableClust = esdTrack->GetTPCNclsF();
 384
 385   Float_t chi2PerCluster = 0.;
 386   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
 387
 388   Float_t clustPerFindClust = 0.;
 389   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
 390
 391   Float_t qpt = 0;
 392   if( fabs(pt)>0 ) qpt = q/fabs(pt);
 393
 394   //
 395   // select primaries
 396   //
 397   Double_t dcaToVertex = -1;
 398   if( fCutsRC->GetDCAToVertex2D() )
 399   {
 400       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 401   }
 402   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 403   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 404   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 405
 406   Double_t vTPCTrackHisto[10] = {nClust,chi2PerCluster,clustPerFindClust,dca[0],dca[1],eta,phi,pt,qpt,vertStatus};
 407   fTPCTrackHisto->Fill(vTPCTrackHisto);
 408
 409   //
 410   // Fill rec vs MC information
 411   //
 412   if(!stack) return;
 413 }
 414
 415
 416 //_____________________________________________________________________________
 417 void AliPerformanceTPC::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/, AliESDEvent* const /*esdEvent*/)
 418 {
 419   // Fill comparison information (constarained parameters)
 420   AliDebug(AliLog::kWarning, "Warning: Not implemented");
 421 }
 422
 423
 424 //_____________________________________________________________________________
 425 void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 426 {
 427   // Process comparison information
 428   //
 429
 430   if(!esdEvent)
 431   {
 432     Error("Exec","esdEvent not available");
 433     return;
 434   }
 435   AliHeader* header = 0;
 436   AliGenEventHeader* genHeader = 0;
 437   AliStack* stack = 0;
 438   TArrayF vtxMC(3);
 439
 440   if(bUseMC)
 441   {
 442     if(!mcEvent) {
 443       Error("Exec","mcEvent not available");
 444       return;
 445     }
 446     // get MC event header
 447     header = mcEvent->Header();
 448     if (!header) {
 449       Error("Exec","Header not available");
 450       return;
 451     }
 452     // MC particle stack
 453     stack = mcEvent->Stack();
 454     if (!stack) {
 455       Error("Exec","Stack not available");
 456       return;
 457     }
 458     // get MC vertex
 459     genHeader = header->GenEventHeader();
 460     if (!genHeader) {
 461       Error("Exec","Could not retrieve genHeader from Header");
 462       return;
 463     }
 464     genHeader->PrimaryVertex(vtxMC);
 465   }
 466
 467   // trigger
 468   if(!bUseMC &&GetTriggerClass()) {
 469     Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
 470     if(!isEventTriggered) return;
 471   }
 472
 473   // get TPC event vertex
 474   const AliESDVertex *vtxESD = NULL;
 475   if(fUseTrackVertex) {
 476     vtxESD = esdEvent->GetPrimaryVertexTracks();
 477   } else {
 478     vtxESD = esdEvent->GetPrimaryVertexTPC();
 479   }
 480   if(!vtxESD) return;
 481
 482   //  events with rec. vertex
 483   Int_t mult=0; Int_t multP=0; Int_t multN=0;
 484
 485   // changed to take all events but store vertex status
 486 //  if(vtxESD->GetStatus() >0)
 487 //  {
 488   // store vertex status
 489   Bool_t vertStatus = vtxESD->GetStatus();
 490   //  Process ESD events
 491   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 492   {
 493     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 494     if(!track) continue;
 495
 496     // if not fUseKinkDaughters don't use tracks with kink index > 0
 497     if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
 498
 499     if(bUseESDfriend && esdFriend && esdFriend->TestSkipBit()==kFALSE)
 500     {
 501       AliESDfriendTrack *friendTrack=esdFriend->GetTrack(iTrack);
 502       if(friendTrack)
 503       {
 504         //
 505         TObject *calibObject=0;
 506         AliTPCseed *seed=0;
 507         for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j) {
 508             if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) {
 509             break;
 510           }
 511         }
 512
 513         //
 514         for (Int_t irow=0;irow<159;irow++) {
 515         if(!seed) continue;
 516
 517           AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 518           if (!cluster) continue;
 519
 520              Float_t gclf[3];
 521              cluster->GetGlobalXYZ(gclf);
 522
 523              //Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 524              //Int_t i[1]={cluster->GetDetector()};
 525              //transform->Transform(x,i,0,1);
 526              //printf("gx %f gy  %f  gz %f \n", cluster->GetX(), cluster->GetY(),cluster->GetZ());
 527              //printf("gclf[0] %f gclf[1]  %f  gclf[2] %f \n", gclf[0], gclf[1],  gclf[2]);
 528
 529              Int_t TPCside;
 530              if(gclf[2]>0.) TPCside=0; // A side
 531              else TPCside=1;
 532
 533              //
 534              //Double_t vTPCClust1[3] = { gclf[0], gclf[1],  TPCside };
 535              //fTPCClustHisto1->Fill(vTPCClust1);
 536
 537              //
 538              Double_t phi = TMath::ATan2(gclf[1],gclf[0]);
 539              if(phi < 0) phi += 2.*TMath::Pi();
 540
 541              //Float_t pad = cluster->GetPad();
 542              //Int_t detector = cluster->GetDetector();
 543              //Double_t vTPCClust[6] = { irow, phi, TPCside, pad, detector, gclf[2] };
 544              Double_t vTPCClust[3] = { irow, phi, TPCside };
 545              fTPCClustHisto->Fill(vTPCClust);
 546         }
 547       }
 548     }
 549
 550     if(GetAnalysisMode() == 0) ProcessTPC(stack,track,esdEvent,vertStatus);
 551     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track,esdEvent,vertStatus);
 552     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track,esdEvent);
 553     else {
 554       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 555       return;
 556     }
 557
 558     // TPC only
 559     if(!fUseHLT){
 560       if(GetAnalysisMode() == 0) {
 561         AliESDtrack *tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
 562         if(!tpcTrack) continue;
 563
 564         // track selection
 565         if( fCutsRC->AcceptTrack(tpcTrack) ) {
 566           mult++;
 567           if(tpcTrack->Charge()>0.) multP++;
 568           if(tpcTrack->Charge()<0.) multN++;
 569         }
 570         if(tpcTrack) delete tpcTrack;
 571       }
 572       else {
 573        // track selection
 574         if( fCutsRC->AcceptTrack(track) ) {
 575           mult++;
 576           if(track->Charge()>0.) multP++;
 577           if(track->Charge()<0.) multN++;
 578         }
 579       }
 580     }
 581     else {
 582       if( fCutsRC->AcceptTrack(track) ) {
 583         //Printf("Still here for HLT");
 584         mult++;
 585         if(track->Charge()>0.) multP++;
 586         if(track->Charge()<0.) multN++;
 587       }
 588     }
 589   }
 590
 591   Double_t vTPCEvent[7] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),mult,multP,multN,vtxESD->GetStatus()};
 592   fTPCEventHisto->Fill(vTPCEvent);
 593 }
 594
 595
 596 //_____________________________________________________________________________
 597 void AliPerformanceTPC::Analyse()
 598 {
 599     //
 600     // Analyse comparison information and store output histograms
 601     // in the folder "folderTPC"
 602     //
 603     TH1::AddDirectory(kFALSE);
 604     TH1::SetDefaultSumw2(kFALSE);
 605     TObjArray *aFolderObj = new TObjArray;
 606     //aFolderObj->SetOwner(); // objects are owned by fanalysisFolder
 607     TString selString;
 608
 609     //
 610     // Cluster histograms
 611     //
 612     AddProjection(aFolderObj, "clust", fTPCClustHisto, 0, 1, 2);
 613
 614     selString = "all";
 615     for(Int_t i=0; i <= 2; i++) {
 616         AddProjection(aFolderObj, "clust", fTPCClustHisto, i, &selString);
 617     }
 618
 619     //fTPCClustHisto->GetAxis(2)->SetRange(1,1); // A-side
 620     //selString = "A_side";
 621     //AddProjection(aFolderObj, fTPCClustHisto, 0, 1, &selString);
 622
 623     //fTPCClustHisto->GetAxis(2)->SetRange(2,2); // C-side
 624     //selString = "C_side";
 625     //AddProjection(aFolderObj, fTPCClustHisto, 0, 1, &selString);
 626
 627     //reset range
 628     fTPCClustHisto->GetAxis(2)->SetRange(1,2);
 629
 630     //
 631     // event histograms
 632     //
 633     for(Int_t i=0; i<=6; i++) {
 634       AddProjection(aFolderObj, "event", fTPCEventHisto, i);
 635     }
 636     AddProjection(aFolderObj, "event", fTPCEventHisto, 4, 5);
 637     AddProjection(aFolderObj, "event", fTPCEventHisto, 0, 1);
 638     AddProjection(aFolderObj, "event", fTPCEventHisto, 0, 3);
 639     AddProjection(aFolderObj, "event", fTPCEventHisto, 1, 3);
 640     AddProjection(aFolderObj, "event", fTPCEventHisto, 2, 3);
 641
 642     // reconstructed vertex status > 0
 643     fTPCEventHisto->GetAxis(6)->SetRange(2,2);
 644     selString = "recVertex";
 645     for(Int_t i=0; i<=5; i++) {
 646       AddProjection(aFolderObj, "event", fTPCEventHisto, i, &selString);
 647     }
 648     AddProjection(aFolderObj, "event", fTPCEventHisto, 4, 5, &selString);
 649     AddProjection(aFolderObj, "event", fTPCEventHisto, 0, 1, &selString);
 650     AddProjection(aFolderObj, "event", fTPCEventHisto, 0, 3, &selString);
 651     AddProjection(aFolderObj, "event", fTPCEventHisto, 1, 3, &selString);
 652     AddProjection(aFolderObj, "event", fTPCEventHisto, 2, 3, &selString);
 653
 654     // reset cuts
 655     fTPCEventHisto->GetAxis(6)->SetRange(1,2);
 656
 657     //
 658     // Track histograms
 659     //
 660     // all with vertex
 661     fTPCTrackHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
 662     fTPCTrackHisto->GetAxis(9)->SetRangeUser(0.5,1.5);
 663     selString = "all_recVertex";
 664     for(Int_t i=0; i <= 9; i++) {
 665         AddProjection(aFolderObj, "track", fTPCTrackHisto, i, &selString);
 666     }
 667
 668      AddProjection(aFolderObj, "track", fTPCTrackHisto, 5, 8, &selString);
 669
 670     for(Int_t i=0; i <= 4; i++) {
 671         AddProjection(aFolderObj, "track", fTPCTrackHisto, i, 5, 7, &selString);
 672     }
 673
 674
 675
 676     // Track histograms (pos with vertex)
 677     fTPCTrackHisto->GetAxis(8)->SetRangeUser(0,1.5);
 678     selString = "pos_recVertex";
 679     for(Int_t i=0; i <= 9; i++) {
 680         AddProjection(aFolderObj, "track", fTPCTrackHisto, i, &selString);
 681     }
 682     for(Int_t i=0; i <= 4; i++) { for(Int_t j=5; j <= 5; j++) { for(Int_t k=j+1; k <= 7; k++) {
 683         AddProjection(aFolderObj, "track", fTPCTrackHisto, i, j, k, &selString);
 684     }  }  }
 685     AddProjection(aFolderObj, "track", fTPCTrackHisto, 0, 1, 2, &selString);
 686     AddProjection(aFolderObj, "track", fTPCTrackHisto, 0, 1, 5, &selString);
 687     AddProjection(aFolderObj, "track", fTPCTrackHisto, 0, 2, 5, &selString);
 688     AddProjection(aFolderObj, "track", fTPCTrackHisto, 1, 2, 5, &selString);
 689     AddProjection(aFolderObj, "track", fTPCTrackHisto, 3, 4, 5, &selString);
 690     AddProjection(aFolderObj, "track", fTPCTrackHisto, 5, 6, 7, &selString);
 691
 692     // Track histograms (neg with vertex)
 693     fTPCTrackHisto->GetAxis(8)->SetRangeUser(-1.5,0);
 694     selString = "neg_recVertex";
 695     for(Int_t i=0; i <= 9; i++) {
 696         AddProjection(aFolderObj, "track", fTPCTrackHisto, i, &selString);
 697     }
 698     for(Int_t i=0; i <= 4; i++) { for(Int_t j=5; j <= 5; j++) { for(Int_t k=j+1; k <= 7; k++) {
 699         AddProjection(aFolderObj, "track", fTPCTrackHisto, i, j, k, &selString);
 700     }  }  }
 701     AddProjection(aFolderObj, "track", fTPCTrackHisto, 0, 1, 2, &selString);
 702     AddProjection(aFolderObj, "track", fTPCTrackHisto, 0, 1, 5, &selString);
 703     AddProjection(aFolderObj, "track", fTPCTrackHisto, 0, 2, 5, &selString);
 704     AddProjection(aFolderObj, "track", fTPCTrackHisto, 1, 2, 5, &selString);
 705     AddProjection(aFolderObj, "track", fTPCTrackHisto, 3, 4, 5, &selString);
 706     AddProjection(aFolderObj, "track", fTPCTrackHisto, 5, 6, 7, &selString);
 707
 708     //restore cuts
 709     fTPCTrackHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
 710     fTPCTrackHisto->GetAxis(9)->SetRangeUser(-0.5,1.5);
 711
 712
 713     printf("exportToFolder\n");
 714     // export objects to analysis folder
 715     fAnalysisFolder = ExportToFolder(aFolderObj);
 716     if (fFolderObj) delete fFolderObj;
 717     fFolderObj = aFolderObj;
 718     aFolderObj=0;
 719 }
 720
 721
 722 //_____________________________________________________________________________
 723 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array)
 724 {
 725   // recreate folder avery time and export objects to new one
 726   //
 727   AliPerformanceTPC * comp=this;
 728   TFolder *folder = comp->GetAnalysisFolder();
 729
 730   TString name, title;
 731   TFolder *newFolder = 0;
 732   Int_t i = 0;
 733   Int_t size = array->GetSize();
 734
 735   if(folder) {
 736      // get name and title from old folder
 737      name = folder->GetName();
 738      title = folder->GetTitle();
 739
 740          // delete old one
 741      delete folder;
 742
 743          // create new one
 744      newFolder = CreateFolder(name.Data(),title.Data());
 745      newFolder->SetOwner();
 746
 747          // add objects to folder
 748      while(i < size) {
 749            newFolder->Add(array->At(i));
 750            i++;
 751          }
 752   }
 753
 754 return newFolder;
 755 }
 756
 757 //_____________________________________________________________________________
 758 Long64_t AliPerformanceTPC::Merge(TCollection* const list)
 759 {
 760   // Merge list of objects (needed by PROOF)
 761
 762   if (!list)
 763   return 0;
 764
 765   if (list->IsEmpty())
 766   return 1;
 767
 768   TIterator* iter = list->MakeIterator();
 769   TObject* obj = 0;
 770   TObjArray* objArrayList = 0;
 771   objArrayList = new TObjArray();
 772
 773   // collection of generated histograms
 774   Int_t count=0;
 775   while((obj = iter->Next()) != 0)
 776   {
 777     AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
 778     if (entry == 0) continue;
 779     if (fgMergeTHnSparse) {
 780         if ((fTPCClustHisto) && (entry->fTPCClustHisto)) { fTPCClustHisto->Add(entry->fTPCClustHisto); }
 781         if ((fTPCEventHisto) && (entry->fTPCEventHisto)) { fTPCEventHisto->Add(entry->fTPCEventHisto); }
 782         if ((fTPCTrackHisto) && (entry->fTPCTrackHisto)) { fTPCTrackHisto->Add(entry->fTPCTrackHisto); }
 783     }
 784     // the analysisfolder is only merged if present
 785     if (entry->fFolderObj) { objArrayList->Add(entry->fFolderObj); }
 786
 787     count++;
 788   }
 789   if (fFolderObj) { fFolderObj->Merge(objArrayList); }
 790   // to signal that track histos were not merged: reset
 791   if (!fgMergeTHnSparse) { fTPCTrackHisto->Reset(); fTPCClustHisto->Reset(); fTPCEventHisto->Reset(); }
 792   // delete
 793   if (objArrayList)  delete objArrayList;  objArrayList=0;
 794 return count;
 795 }
 796
 797
 798 //_____________________________________________________________________________
 799 TFolder* AliPerformanceTPC::CreateFolder(TString name, TString title)
 800 {
 801 // create folder for analysed histograms
 802 //
 803 TFolder *folder = 0;
 804   folder = new TFolder(name.Data(),title.Data());
 805
 806   return folder;
 807 }
 808
 809 //_____________________________________________________________________________
 810 TTree* AliPerformanceTPC::CreateSummary()
 811 {
 812     // implementaion removed, switched back to use AliPerformanceSummary (now called in AliPerformanceTask)
 813     return 0;
 814 }
 815