PWG1/AliComparisonDEdx.cxx

   1 //------------------------------------------------------------------------------\r
   2 // Implementation of AliComparisonDEdx class. It keeps information from \r
   3 // comparison of reconstructed and MC particle tracks. In addtion, \r
   4 // it keeps selection cuts used during comparison. The comparison \r
   5 // information is stored in the ROOT histograms. Analysis of these \r
   6 // histograms can be done by using Analyse() class function. The result of \r
   7 // the analysis (histograms) are stored in the output picture_dedx.root file.\r
   8 //  \r
   9 // Author: J.Otwinowski 04/02/2008 \r
  10 //------------------------------------------------------------------------------\r
  11 \r
  12 \r
  13 /*\r
  14   //after running analysis, read the file, and get component\r
  15   gSystem->Load("libPWG1.so");\r
  16   TFile f("Output.root");\r
  17   AliComparisonDEdx * comp = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");\r
  18 \r
  19   // analyse comparison data (output stored in pictures_dedx.root)\r
  20   comp->Analyse();\r
  21 \r
  22   // TPC track length parameterisation\r
  23   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function\r
  24   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);\r
  25   fl2.SetParameter(1,1);\r
  26   fl2.SetParameter(0,1);\r
  27 */\r
  28 \r
  29 #include <iostream>\r
  30 \r
  31 #include "TFile.h"\r
  32 #include "TCint.h"\r
  33 #include "TH3F.h"\r
  34 #include "TH2F.h"\r
  35 #include "TF1.h"\r
  36 #include "TProfile.h"\r
  37 #include "TProfile2D.h"\r
  38 #include "TGraph2D.h"\r
  39 #include "TCanvas.h"\r
  40 #include "TGraph.h"\r
  41 //\r
  42 #include "AliESDEvent.h"\r
  43 #include "AliESD.h"\r
  44 #include "AliESDfriend.h"\r
  45 #include "AliESDfriendTrack.h"\r
  46 #include "AliRecInfoCuts.h" \r
  47 #include "AliMCInfoCuts.h" \r
  48 #include "AliLog.h" \r
  49 //\r
  50 #include "AliMathBase.h"\r
  51 #include "AliTreeDraw.h"\r
  52 //#include "TStatToolkit.h"\r
  53 \r
  54 #include "AliMCInfo.h" \r
  55 #include "AliESDRecInfo.h" \r
  56 #include "AliComparisonDEdx.h" \r
  57 \r
  58 using namespace std;\r
  59 \r
  60 ClassImp(AliComparisonDEdx)\r
  61 \r
  62 //_____________________________________________________________________________\r
  63 AliComparisonDEdx::AliComparisonDEdx():\r
  64   TNamed("AliComparisonDEdx","AliComparisonDEdx"),\r
  65 \r
  66   // dEdx \r
  67   fTPCSignalNormTan(0), \r
  68   fTPCSignalNormSPhi(0),\r
  69   fTPCSignalNormTPhi(0), \r
  70   //\r
  71   fTPCSignalNormTanSPhi(0),\r
  72   fTPCSignalNormTanTPhi(0),\r
  73   fTPCSignalNormTanSPt(0), \r
  74   \r
  75   // Cuts \r
  76   fCutsRC(0), \r
  77   fCutsMC(0),\r
  78   fMCPtMin(0),\r
  79   fMCAbsTanThetaMax(0),\r
  80   fMCPdgCode(0)\r
  81 {\r
  82   InitHisto();\r
  83   InitCuts();\r
  84 }\r
  85 \r
  86 //_____________________________________________________________________________\r
  87 AliComparisonDEdx::~AliComparisonDEdx(){\r
  88    \r
  89   if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0; \r
  90   if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;\r
  91   if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;\r
  92   //\r
  93   if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;\r
  94   if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;\r
  95   if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;\r
  96 }\r
  97 \r
  98 //_____________________________________________________________________________\r
  99 void AliComparisonDEdx::InitHisto()\r
 100 {\r
 101   // Init histograms\r
 102   \r
 103   // TPC dEdx\r
 104   fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70); \r
 105   fTPCSignalNormTan->SetXTitle("tan(#theta)");\r
 106   fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");\r
 107 \r
 108   fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);\r
 109   fTPCSignalNormSPhi->SetXTitle("sin(#phi)");\r
 110   fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");\r
 111 \r
 112   fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); \r
 113   fTPCSignalNormTPhi->SetXTitle("tan(#phi)");\r
 114   fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");\r
 115 \r
 116   fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);\r
 117   fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");\r
 118   fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");\r
 119   fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");\r
 120 \r
 121   fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);\r
 122   fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");\r
 123   fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");\r
 124   fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");\r
 125 \r
 126   fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); \r
 127   fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");\r
 128   fTPCSignalNormTanSPt->SetYTitle("#sqrt{p_{t}}");\r
 129   fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");\r
 130 }\r
 131 \r
 132 void AliComparisonDEdx::InitCuts()\r
 133 {\r
 134   // Init cuts\r
 135   if(!fCutsMC) \r
 136     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");\r
 137   if(!fCutsRC) \r
 138     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");\r
 139 }\r
 140 \r
 141 //_____________________________________________________________________________\r
 142 void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){\r
 143 \r
 144   // Fill dE/dx  comparison information\r
 145   \r
 146   Float_t mcpt = infoMC->GetParticle().Pt();\r
 147   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);\r
 148   Float_t mprim = infoMC->GetPrim();\r
 149 \r
 150   // distance to Prim. vertex \r
 151   const Double_t* dv = infoMC->GetVDist(); \r
 152 \r
 153   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();\r
 154   \r
 155   // Check selection cuts \r
 156   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; \r
 157   if (!isPrim) return;\r
 158   if (infoRC->GetStatus(1)!=3) return;\r
 159   if (!infoRC->GetESDtrack()) return;  \r
 160   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;\r
 161   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;\r
 162   //if (mprim>1.4) return;\r
 163   //if (mprim<0.5) return;\r
 164   if (mprim > fCutsMC->GetMaxTPCSignal()) return;\r
 165   if (mprim < fCutsMC->GetMinTPCSignal()) return;\r
 166   if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;\r
 167   //\r
 168   Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();\r
 169   Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();\r
 170   Float_t tphi =  sphi/TMath::Sqrt(1-sphi*sphi);\r
 171 \r
 172   if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;\r
 173   //if (mcpt>0.5) {\r
 174   if (mcpt > GetMCPtMin()) {\r
 175     fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset\r
 176   }\r
 177 \r
 178   //if (TMath::Abs(tantheta)<0.5){\r
 179   if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){\r
 180     fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset\r
 181     fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset\r
 182   }\r
 183   fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);    \r
 184   fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);    \r
 185   fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);\r
 186 }\r
 187 \r
 188 //_____________________________________________________________________________\r
 189 Long64_t AliComparisonDEdx::Merge(TCollection* list) \r
 190 {\r
 191   // Merge list of objects (needed by PROOF)\r
 192 \r
 193   if (!list)\r
 194   return 0;\r
 195 \r
 196   if (list->IsEmpty())\r
 197   return 1;\r
 198 \r
 199   TIterator* iter = list->MakeIterator();\r
 200   TObject* obj = 0;\r
 201 \r
 202   // collection of generated histograms\r
 203   Int_t count=0;\r
 204   while((obj = iter->Next()) != 0) \r
 205   {\r
 206     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);\r
 207     if (entry == 0) { \r
 208       Error("Add","Attempt to add object of class: %s to a %s",\r
 209           obj->ClassName(),this->ClassName());\r
 210       return -1;\r
 211   }\r
 212 \r
 213   fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);\r
 214   fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);\r
 215   fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);\r
 216   //\r
 217   fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);\r
 218   fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);\r
 219   fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);\r
 220 \r
 221   count++;\r
 222   }\r
 223 \r
 224 return count;\r
 225 }\r
 226 \r
 227 //_____________________________________________________________________________\r
 228 void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)\r
 229 {\r
 230   // Process comparison information\r
 231   Process(infoMC,infoRC);\r
 232 }\r
 233 \r
 234 //_____________________________________________________________________________\r
 235 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)\r
 236 {\r
 237   // Make resolution histograms\r
 238   TH1F *hisr, *hism;\r
 239   if (!gPad) new TCanvas;\r
 240   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);\r
 241   if (type) return hism;\r
 242   else \r
 243     return hisr;\r
 244 }\r
 245 \r
 246 //_____________________________________________________________________________\r
 247 void AliComparisonDEdx::Analyse()\r
 248 {\r
 249   // Analyse output histograms\r
 250   \r
 251   AliComparisonDEdx * comp=this;\r
 252   TH1F *hiss=0;\r
 253   TGraph2D * gr=0;\r
 254 \r
 255   TFile *fp = new TFile("pictures_dedx.root","recreate");\r
 256   fp->cd();\r
 257 \r
 258   TCanvas * c = new TCanvas("TPCdedx","TPC dedx");\r
 259   c->cd();\r
 260 \r
 261   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);\r
 262   hiss->SetXTitle("Tan(#theta)");\r
 263   hiss->SetYTitle("#sigma_{dEdx}");\r
 264   hiss->Draw();\r
 265   hiss->Write("TPCdEdxResolTan");\r
 266   //\r
 267   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1); \r
 268   hiss->SetXTitle("Tan(#theta)");\r
 269   hiss->SetYTitle("<dEdx>");\r
 270   hiss->Draw(); \r
 271   hiss->Write("TPCdEdxMeanTan");\r
 272   //\r
 273   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);\r
 274   gr->GetXaxis()->SetTitle("Tan(#theta)");\r
 275   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");\r
 276   gr->GetZaxis()->SetTitle("<dEdx>");\r
 277   gr->Draw("colz"); \r
 278   gr->GetHistogram()->Write("TPCdEdxMeanTanPt_1");\r
 279   //\r
 280   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);\r
 281   gr->GetXaxis()->SetTitle("Tan(#theta)");\r
 282   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");\r
 283   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");\r
 284   gr->Draw("colz"); \r
 285   gr->GetHistogram()->Write("TPCdEdxMeanTanPt_2");\r
 286 \r
 287   fp->Close();\r
 288 }\r