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