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