TRD/Macros/AliTRDhistTRDpid.C

   1 // ////////////////////////////////////////////////////////////////////////////
   2 //   Macro to get TRD signal distribution from all six layers of TRD
   3 //   This function takes the data from following directories
   4 //   "rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5",//   "rmom6", "rmom8", "rmom10"
   5 //    corresponding to momenta
   6 //   0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0
   7 //   Running this macro will generate the root file TRDdEdxHistogramsV1.root
   8 //   containing histograms for all momenta.
   9 //   Prashant Shukla (shukla@physi.uni-heidelberg.de)
  10 // ////////////////////////////////////////////////////////////////////////////
  11
  12 #if !defined( __CINT__) || defined(__MAKECINT__)
  13 #include <TFile.h>
  14 #include <TH1.h>
  15 #include <TH2.h>
  16 #include <TCanvas.h>
  17 #include <TVector3.h>
  18 #include <TPDGCode.h>
  19 #include <TParticle.h>
  20 #include <TTree.h>
  21
  22 #include "AliTRDtrack.h"
  23 #include "AliRunLoader.h"
  24 #include "AliLoader.h"
  25 #include "AliESD.h"
  26 #include "AliRun.h"
  27 #include "AliStack.h"
  28 #include "AliHeader.h"
  29 #include "AliGenEventHeader.h"
  30 #else
  31 #endif
  32
  33
  34 //_____________________________________________________________________________
  35
  36 Bool_t AliTRDhistTRDpid(Int_t oFile=1)
  37 {
  38   //
  39   //
  40
  41   const Int_t kNMom=11;
  42
  43   TH1F *h1dEdxEL[kNMom];
  44   TH1F *h1dEdxPI[kNMom];
  45   TH1F *h1dEdxMU[kNMom];
  46   TH1F *h1dEdxKA[kNMom];
  47   TH1F *h1dEdxPR[kNMom];
  48   TH1F *h1MaxTimBinEL[kNMom];
  49   TH1F *h1MaxTimBinPI[kNMom];
  50   Char_t text[200];
  51
  52   // Histograms for dedx of e and pi track in TRD
  53
  54   for (Int_t imom = 0; imom < kNMom; imom++) {
  55     sprintf(text,"h1dEdxEL%01d",imom+1);
  56     h1dEdxEL[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
  57     h1dEdxEL[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
  58     h1dEdxEL[imom]->GetYaxis()->SetTitle("Entries");
  59     h1dEdxEL[imom]->SetLineColor(kRed);
  60
  61     sprintf(text,"h1dEdxPI%01d",imom+1);
  62     h1dEdxPI[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
  63     //  TH1F *h1dEdxPI[imom] = new TH1F("","",250,0,2500);
  64     h1dEdxPI[imom]->GetXaxis()->SetTitle("dE/dx_{TRD} (arbitrary units)");
  65     h1dEdxPI[imom]->GetYaxis()->SetTitle("Entries");
  66     h1dEdxPI[imom]->SetLineColor(kBlue);
  67
  68     sprintf(text,"h1dEdxMU%01d",imom+1);
  69     h1dEdxMU[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
  70     h1dEdxMU[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
  71     h1dEdxMU[imom]->GetYaxis()->SetTitle("Entries");
  72     h1dEdxMU[imom]->SetLineColor(kGreen);
  73
  74     sprintf(text,"h1dEdxKA%01d",imom+1);
  75     h1dEdxKA[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
  76     h1dEdxKA[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
  77     h1dEdxKA[imom]->GetYaxis()->SetTitle("Entries");
  78     h1dEdxKA[imom]->SetLineColor(7);
  79
  80     sprintf(text,"h1dEdxPR%01d",imom+1);
  81     h1dEdxPR[imom] = new TH1F(text,"dE/dx  Dist.",200,0,4000);
  82     h1dEdxPR[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)");
  83     h1dEdxPR[imom]->GetYaxis()->SetTitle("Entries");
  84     h1dEdxPR[imom]->SetLineColor(6);
  85
  86     sprintf(text,"h1MaxTimBinEL%01d",imom+1);
  87     h1MaxTimBinEL[imom] = new TH1F(text,"Dist. Time Bin of max. Cluster for e(Red) and pi(Blue)",30,0,30);
  88     h1MaxTimBinEL[imom]->GetXaxis()->SetTitle("Time Bin of Maximum Cluster");
  89     h1MaxTimBinEL[imom]->GetYaxis()->SetTitle("Entries");
  90     h1MaxTimBinEL[imom]->SetLineColor(2);
  91
  92     sprintf(text,"h1MaxTimBinPI%01d",imom+1);
  93     h1MaxTimBinPI[imom] = new TH1F(text,"Time Bin of max. Cluster Pion",30,0,30);
  94     h1MaxTimBinPI[imom]->SetLineColor(4);
  95   }
  96
  97   // PID
  98   Int_t partCode[5] =
  99     {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
 100
 101   // Files
 102
 103   Char_t *FileDir[11] = {"rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5", "rmom6", "rmom8", "rmom10"};
 104
 105   Float_t mome[11]={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0};
 106
 107   Char_t gAliceFileName[200];
 108   Char_t esdFileName[200];
 109   for (Int_t imom = 0; imom < kNMom; imom++) {
 110
 111     sprintf(gAliceFileName,"~/work/data/%s/galice.root", FileDir[imom]);
 112     sprintf(esdFileName,"~/work/data/%s/AliESDs.root", FileDir[imom]);
 113
 114   // open run loader and load gAlice, kinematics and header
 115   AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName);
 116   if (!runLoader) {
 117     Error("CheckESD", "getting run loader from file %s failed",gAliceFileName);
 118     return kFALSE;
 119   }
 120   runLoader->LoadgAlice();
 121   gAlice = runLoader->GetAliRun();
 122   if (!gAlice) {
 123     Error("CheckESD", "no galice object found");
 124     return kFALSE;
 125   }
 126   runLoader->LoadKinematics();
 127   runLoader->LoadHeader();
 128
 129   //  AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
 130   // open the ESD file
 131   TFile* esdFile = TFile::Open(esdFileName);
 132   if (!esdFile || !esdFile->IsOpen()) {
 133     Error("CheckESD", "opening ESD file %s failed",esdFileName);
 134     return kFALSE;
 135   }
 136   AliESD* esd = new AliESD;
 137   TTree* tree = (TTree*) esdFile->Get("esdTree");
 138   if (!tree) {
 139     Error("CheckESD", "no ESD tree found");
 140     return kFALSE;
 141   }
 142   tree->SetBranchAddress("ESD", &esd);
 143
 144   // loop over events
 145   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
 146     runLoader->GetEvent(iEvent);
 147     AliStack* stack = gAlice->Stack();
 148     Int_t nParticles = stack->GetNtrack();
 149     TArrayF vertex(3);
 150     runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 151     Int_t nGene=0, nGenpi=0;
 152     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 153       TParticle* particle = stack->Particle(iParticle);
 154       if (!particle) continue;
 155       if (particle->Pt() < 0.001) continue;
 156       //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
 157       TVector3 dVertex(particle->Vx() - vertex[0],
 158                        particle->Vy() - vertex[1],
 159                        particle->Vz() - vertex[2]);
 160       if (dVertex.Mag() > 0.0001) continue;
 161       if (TMath::Abs(particle->GetPdgCode()) == partCode[0]) nGene++;
 162       if (TMath::Abs(particle->GetPdgCode()) == partCode[2]) nGenpi++;
 163     }  // end loop over particles
 164     printf("Gen Prim. el= %d, Gen Prim. pi= %d in Event No.= %d in File = %s \n", nGene, nGenpi, iEvent, FileDir[imom]);
 165
 166     // get the event summary data
 167     tree->GetEvent(iEvent);
 168     if (!esd) {
 169       Error("CheckESD", "no ESD object found for event %d", iEvent);
 170       return kFALSE;
 171     }
 172
 173     ///////////////////////////////////////////
 174     //    AliTRDpidESD::MakePID(esd);
 175     //    AliESDpid::MakePID(esd);
 176     ///////////////////////////////////////////
 177
 178     // Initializing number of electrons and pions
 179     Int_t nne=0, nnpi=0;
 180     // loop over tracks
 181     Int_t nTracks = esd->GetNumberOfTracks();
 182     printf("Found %d tracks. \n",nTracks);
 183
 184     for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 185       AliESDtrack* track = esd->GetTrack(iTrack);
 186       // select tracks of selected particles
 187       Int_t label = TMath::Abs(track->GetLabel());
 188       if (label > stack->GetNtrack()) continue;     // background
 189       TParticle* particle = stack->Particle(label);
 190       //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
 191       TVector3 dVertex(particle->Vx() - vertex[0],
 192                        particle->Vy() - vertex[1],
 193                        particle->Vz() - vertex[2]);
 194       if (dVertex.Mag() > 0.0001) continue;
 195       if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
 196       if (track->GetConstrainedChi2() > 1e9) continue;
 197       Double_t p[3];
 198       track->GetConstrainedPxPyPz(p);
 199       TVector3 pTrack(p);
 200       //      Double_t mom = track->GetP();
 201
 202       // PID/////////////////////////////////////////////////////
 203       if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
 204       Int_t iGen = 5;
 205       for (Int_t i = 0; i < 5; i++) {
 206         if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;
 207       }
 208
 209       // Filling muon kan and proton//////////////////////////////
 210       Int_t Itot=0;
 211       if(track->GetTRDsignal()) Itot=1;
 212       if(!Itot) continue;
 213
 214       // dE/dx and TEST TRD
 215       Double_t TRDsignal = track->GetTRDsignal();
 216       //      if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal);
 217       //      if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal);
 218       //      if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal);
 219       //      if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal);
 220       //      if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal);
 221
 222       //Filling elcron and pions////////////////////////////////
 223       const Int_t kNPlane = 6;
 224       // Track Selection
 225       Double_t Ise=1.;
 226       for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
 227         Ise*=track->GetTRDsignals(iPlane);
 228       }
 229       if(Ise==0) continue;
 230
 231       for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
 232         Double_t TRDsignal = track->GetTRDsignals(iPlane);
 233         if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal);
 234         if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal);
 235         if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal);
 236         if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal);
 237         if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal);
 238
 239         if (iGen == 0) h1MaxTimBinEL[imom]->Fill(track->GetTRDTimBin(iPlane));
 240         if (iGen == 2) h1MaxTimBinPI[imom]->Fill(track->GetTRDTimBin(iPlane));
 241       }
 242
 243       //      Int_t TRDncls = track->GetTRDncls();
 244       if (iGen == 0) nne++;
 245       if (iGen == 2) nnpi++;
 246     } // end of loop over tracks
 247     printf("Electron Tracks = %d, Pion Tracks =  %d \n", nne, nnpi);
 248   }// end of loop over events
 249   delete esd;
 250   esdFile->Close();
 251   delete esdFile;
 252
 253   runLoader->UnloadHeader();
 254   runLoader->UnloadKinematics();
 255   delete runLoader;
 256   } // Loop over files
 257
 258
 259   Char_t fileprint[200];
 260
 261   // draw the histograms if not in batch mode
 262
 263   new TCanvas;
 264   gPad->SetBorderMode(0);
 265   gPad->SetFillColor(0);
 266   //  gStyle->SetStatColor(0);
 267   h1dEdxPI[0]->DrawCopy();
 268   h1dEdxEL[0]->DrawCopy("SAME");
 269   //  h1dEdxMU[0]->DrawCopy("SAME");
 270   //  h1dEdxKA[0]->DrawCopy("SAME");
 271   //  h1dEdxPR[0]->DrawCopy("SAME");
 272
 273   new TCanvas;
 274   gPad->SetBorderMode(0);
 275   gPad->SetFillColor(0);
 276   h1MaxTimBinEL[0]->DrawCopy();
 277   h1MaxTimBinPI[0]->DrawCopy("SAME");
 278
 279
 280   // write the output histograms to a file
 281
 282   Char_t fileresponse[200];
 283   sprintf(fileresponse,"TRDdEdxHistogramsV%d.root",oFile);
 284
 285   TFile* outputFile = TFile::Open(fileresponse, "recreate");
 286   if (!outputFile || !outputFile->IsOpen()) {
 287     Error("CheckESD", "opening output file check.root failed");
 288     return kFALSE;
 289   }
 290
 291   for (Int_t imom = 0; imom < kNMom; imom++) {
 292     h1dEdxPI[imom]->Write();
 293     h1dEdxEL[imom]->Write();
 294     h1dEdxMU[imom]->Write();
 295     h1dEdxKA[imom]->Write();
 296     h1dEdxPR[imom]->Write();
 297     h1MaxTimBinEL[imom]->Write();
 298     h1MaxTimBinPI[imom]->Write();
 299   }
 300
 301   delete [] h1dEdxPI;
 302   delete [] h1dEdxEL;
 303   delete [] h1dEdxMU;
 304   delete [] h1dEdxKA;
 305   delete [] h1dEdxPR;
 306   delete [] h1MaxTimBinEL;
 307   delete [] h1MaxTimBinPI;
 308   delete outputFile;
 309
 310   // result of check
 311   Info("GetPIDsignals", "GetPIDsignals was successfull");
 312   return kTRUE;
 313
 314 } // end of main program
 315