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