| 1 |
// sample to plot ND data |
| 2 |
// 04/26/2007 |
| 3 |
// |
| 4 |
// |
| 5 |
// C/C++ headers |
| 6 |
// |
| 7 |
#include "TrAng.h" |
| 8 |
#include "QtoInc.h" |
| 9 |
|
| 10 |
#include <fstream> |
| 11 |
#include <iostream> |
| 12 |
#include <stdio.h> |
| 13 |
// |
| 14 |
// ROOT headers |
| 15 |
// |
| 16 |
#include <TTree.h> |
| 17 |
#include <TObject.h> |
| 18 |
#include <TList.h> |
| 19 |
#include <TArrayI.h> |
| 20 |
#include <TSystem.h> |
| 21 |
#include <TSystemDirectory.h> |
| 22 |
#include <TString.h> |
| 23 |
#include <TFile.h> |
| 24 |
#include <TCanvas.h> |
| 25 |
//#include <TMatrixD.h> |
| 26 |
#include "TStyle.h" |
| 27 |
#include <TH2F.h> |
| 28 |
|
| 29 |
#include <PamLevel2.h> |
| 30 |
#include <TStyle.h> |
| 31 |
|
| 32 |
using namespace std; |
| 33 |
|
| 34 |
int main(int argc, char* argv[]){ |
| 35 |
|
| 36 |
string outdir = "/home/pamelaprod/malakhov/QtoInc/Picture/"; |
| 37 |
Bool_t AnglHist = false; |
| 38 |
Bool_t DiffHist = false; |
| 39 |
Bool_t PamEff = false; |
| 40 |
Bool_t DoTr = false; |
| 41 |
Bool_t PamAngTime = false; |
| 42 |
Bool_t PamExp = false; |
| 43 |
|
| 44 |
if(argc<2){ |
| 45 |
cout<<"You have to insert at least a file to analisys \nUsing --help for more information\n"; |
| 46 |
exit(1); |
| 47 |
} |
| 48 |
|
| 49 |
if(!strcmp(argv[1],"--help")){ |
| 50 |
cout<<"\nUsage \n"; |
| 51 |
cout<<"\n QtoInc --WorkDir directory_with_level2_files --OutPuth directory_for_output files [Options] \n"; |
| 52 |
cout<<"\n --WorkDir full puth to the Level2 file \n"; |
| 53 |
cout<<" --OutPuth full puth for putputh files \n"; |
| 54 |
cout<<"\n options are:\n"; |
| 55 |
cout<<"--help print this help and exit \n"; |
| 56 |
cout<<"-DoTr Use DoTrack2 procedure for calculating direction of particle flight into Pamela \n"; |
| 57 |
cout<<"-PamEff Calculating Pamela's angular efficiency \n"; |
| 58 |
cout<<"-DiffHist Getting hystogramm for diference between direction getting using DoTrack2 and using axv & ayv values \n"; |
| 59 |
cout<<"-AngHist Getting hystogramm for various angles(Pitch Angle, Pamela's Main axis angles, etc) [default] \n"; |
| 60 |
cout<<"-PamAngTime Not Useful now \n"; |
| 61 |
cout<<"-PamExp Calculating exposition for each possible direction in Pamela's aperture \n"; |
| 62 |
cout<<"\nFor Example: \n"; |
| 63 |
cout<<"\n./QtoInc.exe --WorkDir /data01/_3/704_3/ --OutPuth /home/pamelaprod/malakhov/QtoInc/Picture5/ -DiffHist -DoTr \n\n"; |
| 64 |
exit(1); |
| 65 |
} |
| 66 |
|
| 67 |
if(!strcmp(argv[3],"--OutPuth")){ |
| 68 |
if (4 >= argc){ |
| 69 |
cout<<"--OutPath needs argument\n See --help\n"; |
| 70 |
}else outdir = argv[4]; |
| 71 |
} |
| 72 |
|
| 73 |
for(Int_t i=4; i<argc; i++){ |
| 74 |
if(!strcmp(argv[i],"-DoTr")) DoTr = true; |
| 75 |
if(!strcmp(argv[i],"-PamEff")) PamEff = true; |
| 76 |
if(!strcmp(argv[i],"-DiffHist")){DiffHist = true; DoTr = true;} |
| 77 |
if(!strcmp(argv[i],"-AngHist")) AnglHist = true; |
| 78 |
if(!strcmp(argv[i],"-PamAngTime")) PamAngTime = true; |
| 79 |
if(!strcmp(argv[i],"-PamExp")) PamExp = true; |
| 80 |
} |
| 81 |
|
| 82 |
if(!PamEff && !DiffHist && !AnglHist && !PamAngTime && !PamExp) AnglHist = true; |
| 83 |
|
| 84 |
if(!strcmp(argv[1],"--WorkDir")){ |
| 85 |
if (2 >= argc){ |
| 86 |
cout<<"--WorkDir needs argument\n See --help\n"; |
| 87 |
exit(1); |
| 88 |
}else Calculate(argv[2],outdir,AnglHist,DiffHist,PamEff,DoTr,PamAngTime,PamExp); |
| 89 |
} |
| 90 |
|
| 91 |
} |
| 92 |
|
| 93 |
//Function to Convert Rigidity to Energy |
| 94 |
|
| 95 |
Float_t ConvertR2T(Float_t R, Float_t M, Float_t Z) |
| 96 |
{ |
| 97 |
// |
| 98 |
// Convert rigidity (GV) in kinetic energy (GeV) per nucleon |
| 99 |
// input = rigidity (GV), mass (Gev/c^2), A,Z. |
| 100 |
// |
| 101 |
|
| 102 |
Float_t EcinperN = (sqrt(pow((Z*R),2)+pow(M,2))-M); |
| 103 |
return EcinperN; |
| 104 |
} |
| 105 |
|
| 106 |
//Function to convert Energy to Rigidity |
| 107 |
|
| 108 |
Float_t ConvertT2R(Float_t T, Float_t M, Float_t A, Float_t Z) |
| 109 |
{ |
| 110 |
// |
| 111 |
// Convert kinetic energy (GeV) per nucleon in rigidity (GV) |
| 112 |
// output = rigidity (GV),input kin energy mass (Gev/c^2), A,Z. |
| 113 |
// |
| 114 |
|
| 115 |
Float_t R= (1/Z)*(sqrt(pow((A*T+M),2)-pow(M,2))); |
| 116 |
return R; |
| 117 |
} |
| 118 |
|
| 119 |
void Calculate(char* dirname, string outdir, Bool_t AnglHist, Bool_t DiffHist, Bool_t PamEff, Bool_t DoTr, Bool_t PamAngTime, Bool_t PamExp){ |
| 120 |
|
| 121 |
/**********************************************/ |
| 122 |
//First initialization, general for all purposes |
| 123 |
/**********************************************/ |
| 124 |
|
| 125 |
//gStyle->SetOptStat(111111); |
| 126 |
TSystemDirectory *workdir = new TSystemDirectory("workdir",dirname); |
| 127 |
TList *flist=workdir->GetListOfFiles(); |
| 128 |
PamLevel2* pam_events = new PamLevel2(); |
| 129 |
PamelaOrientation* PO = new PamelaOrientation(); |
| 130 |
TTree *T = pam_events->GetPamTree(flist,"treename"); |
| 131 |
ULong_t nevents = T->GetEntries(); |
| 132 |
|
| 133 |
cout<<"Number of events: "<<nevents<<endl; |
| 134 |
/********************************************************************************/ |
| 135 |
/*****NEED TO CHANGE FOR OTHER COMPUTERS****************************************/ |
| 136 |
pam_events->GetTrkLevel2()->LoadField("/data01/pamhome/installed/pamela_software/calib/trk-param/field_param-0/"); |
| 137 |
/********************************************************************************/ |
| 138 |
|
| 139 |
Int_t nz = 6; Float_t zin[6]; |
| 140 |
for (Int_t ip=0;ip<nz;ip++) {zin[ip] = pam_events->GetToFLevel2()->GetZTOF(pam_events->GetToFLevel2()->GetToFPlaneID(ip));cout<<zin[ip]<<endl;} |
| 141 |
Trajectory *tr = new Trajectory(nz,zin); |
| 142 |
|
| 143 |
PamTrack *track; |
| 144 |
|
| 145 |
Int_t ntr; |
| 146 |
Float_t Argv = 0; |
| 147 |
Double_t PamAzim = 0; |
| 148 |
Double_t PamZenit = 0; |
| 149 |
Double_t MyAzim = 0; |
| 150 |
Double_t MyZenit = 0; |
| 151 |
Double_t Px; |
| 152 |
Double_t Py; |
| 153 |
Double_t Pz; |
| 154 |
|
| 155 |
/*********************************************************************************/ |
| 156 |
// Angular exposure. How much time Pamela have any Azimutal and Zenital angle |
| 157 |
/*********************************************************************************/ |
| 158 |
|
| 159 |
TH2F *hhigh; |
| 160 |
TH2F *hlow; |
| 161 |
|
| 162 |
TFile fhigh((TString)outdir+"PamAngEfficiencyHighEnergy.root"); |
| 163 |
TFile flow((TString)outdir+"PamAngEfficiencyLowEnergy.root"); |
| 164 |
if(PamExp){ |
| 165 |
|
| 166 |
if (fhigh.IsZombie()||flow.IsZombie()){ |
| 167 |
cout<<"Problem with Hystogrammfiles"<<endl; |
| 168 |
exit(1); |
| 169 |
}else{ |
| 170 |
hhigh = (TH2F*)fhigh.Get("PamAngEffhigh"); |
| 171 |
hlow = (TH2F*)flow.Get("PamAngEfflow"); |
| 172 |
} |
| 173 |
|
| 174 |
} |
| 175 |
|
| 176 |
/*********************************************************************************/ |
| 177 |
//Histogramms for Count and exposition of Angles (Pitch, Pamela's Main axis, etc ) |
| 178 |
/*********************************************************************************/ |
| 179 |
|
| 180 |
//if(AnglHist){ |
| 181 |
|
| 182 |
TH1F *PitchExpositionWithoutBrazil = new TH1F("PitchExpositionWithoutBrazil", "Pitch Exposition without Brazil with Track", 360, 0, 180); |
| 183 |
TH1F *PitchExpositionEquator = new TH1F("PitchExpositionEquator", "Pitch Exposition in Equator with Track", 360, 0, 180); |
| 184 |
TH1F *PitchExpositionBrazil = new TH1F("PitchExpositionBrazil", "Pitch Exposition Brazil with Track", 360, 0, 180); |
| 185 |
TH1F *MainAxisPamelaExpositionWithoutBrazil = new TH1F("MainAxisPamelaPitchExpositionWithoutBrazil", "Main Axis of Pamela Exposition without Brazil with Track", 360, 0, 180); |
| 186 |
TH1F *MainAxisPamelaExpositionEquator = new TH1F("MainAxisPamelaPitchExpositionEquator", "Main Axis of Pamela Exposition in Equator with Track", 360, 0, 180); |
| 187 |
TH1F *MainAxisPamelaExpositionBrazil = new TH1F("MainAxisPamelaExpositionBrazil", "Main Axis of Pamela Exposition Brazil with Track", 360, 0, 180); |
| 188 |
TH1F *PitchExpositionNoOrientationWithoutBrazil = new TH1F("PitchExpositionNoOrientationWithoutBrazil", "Pitch Exposition without orientation Information without Brazil with Track", 360, 0, 180); |
| 189 |
TH1F *PitchExpositionNoOrientationEquator = new TH1F("PitchExpositionNoOrientationEquator", "Pitch Exposition without orientation Information in Equator with Track", 360, 0, 180); |
| 190 |
TH1F *PitchExpositionNoOrientationBrazil = new TH1F("PitchExpositionNoOrientationBrazil", "Pitch Exposition without orientation Information in Brazil with Track", 360, 0, 180); |
| 191 |
TH1F *DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil", "Diference between no orientation data and with orientation without Brazil with Track", 360, 0, 180); |
| 192 |
TH1F *DiferenceBetweenNoOrientationAndWithOrientationEquator = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationEquator", "Diference between no orientation data and with orientation in Equator with Track", 360, 0, 180); |
| 193 |
TH1F *DiferenceBetweenNoOrientationAndWithOrientationBrazil = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationBrazil", "Diference between no orientation data and with orientation in Brazil with Track", 360, 0, 180); |
| 194 |
TH1F *DifferenceCountWithoutBrazil = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3); |
| 195 |
TH1F *DifferenceCountEquator = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3); |
| 196 |
TH1F *DifferenceCountInBrazil = new TH1F("DiferenceCountInBrazil","Diference in Brazil",150,0.001,3); |
| 197 |
|
| 198 |
TH1F *PitchCountWithoutBrazil = new TH1F("PitchCountWithoutBrazil", "Pitch Count without Brazil with Track", 360, 0, 180); |
| 199 |
TH1F *PitchCountEquator = new TH1F("PitchCountEquator", "Pitch Count in Equator with Track", 360, 0, 180); |
| 200 |
TH1F *PitchCountBrazil = new TH1F("PitchCountBrazil", "Pitch Count Brazil with Track", 360, 0, 180); |
| 201 |
TH1F *MainAxisPamelaCountWithoutBrazil = new TH1F("MainAxisPamelaPitchCountWithoutBrazil", "Main Axis of Pamela Count without Brazil with Track", 360, 0, 180); |
| 202 |
TH1F *MainAxisPamelaCountEquator = new TH1F("MainAxisPamelaPitchCountEquator", "Main Axis of Pamela Count in Equator with Track", 360, 0, 180); |
| 203 |
TH1F *MainAxisPamelaCountBrazil = new TH1F("MainAxisPamelaCountBrazil", "Main Axis of Pamela Count Brazil with Track", 360, 0, 180); |
| 204 |
TH1F *PitchCountNoOrientationWithoutBrazil = new TH1F("PitchCountNoOrientationWithoutBrazil", "Pitch Count without orientation Information without Brazil with Track", 360, 0, 180); |
| 205 |
TH1F *PitchCountNoOrientationEquator = new TH1F("PitchCountNoOrientationEquator", "Pitch Count without orientation Information in Equator with Track", 360, 0, 180); |
| 206 |
TH1F *PitchCountNoOrientationBrazil = new TH1F("PitchCountNoOrientationBrazil", "Pitch Count without orientation Information in Brazil with Track", 360, 0, 180); |
| 207 |
//TH1F *Diference2BetweenNoOrientationAndWithOrientationWithoutBrazil = new TH1F("Diference2BetweenNoOrientationAndWithOrientationWithoutBrazil", "Diference between no orientation data and with orientation without Brazil with Track", 360, 0, 180); |
| 208 |
//TH1F *Diference2BetweenNoOrientationAndWithOrientationEquator = new TH1F("Diference2BetweenNoOrientationAndWithOrientationEquator", "Diference between no orientation data and with orientation in Equator with Track", 360, 0, 180); |
| 209 |
//TH1F *Diference2BetweenNoOrientationAndWithOrientationBrazil = new TH1F("Diference2BetweenNoOrientationAndWithOrientationBrazil", "Diference between no orientation data and with orientation in Brazil with Track", 360, 0, 180); |
| 210 |
//TH1F *Difference2CountWithoutBrazil = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3); |
| 211 |
//TH1F *Difference2CountEquator = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3); |
| 212 |
//TH1F *Difference2CountInBrazil = new TH1F("DiferenceCountInBrazil","Diference in Brazil",150,0.001,3); |
| 213 |
//TH1F *EastPitchAllCount = new TH1F("EastPitchAllCount","East Count",) |
| 214 |
|
| 215 |
TH1F *ZenitAngle = new TH1F("ZenitAngle","Zenit Angle",360,120,180); |
| 216 |
//TH1F *ZenitAngleEquator = new TH1F("ZenitAngleWithoutEquator","Zenit Angle in Equator",360,120,180); |
| 217 |
//TH1F *ZenitAngleBrazil = new TH1F("ZenitAngleBrazil","Zenit Angle in Brazil",360,120,180); |
| 218 |
|
| 219 |
TH2F *PitchvsLatInBrazil = new TH2F("PitchvsLatInBrazil","Pitch Angle, deg; Latitude, deg",360,-40,-5,360,40,110); |
| 220 |
TH2F *LieveTimevsLatInBrazil = new TH2F("LieveTimevsLatInBrazil","Live time, deg; Latitude deg",360,-40,-5,1000,0,80); |
| 221 |
|
| 222 |
|
| 223 |
TCanvas *c1 = new TCanvas("c1","Exposition and count in Brazil",1200,3000); |
| 224 |
TCanvas *c2 = new TCanvas("c2","Exposition and Count without Brazil",1200,3000); |
| 225 |
TCanvas *c3 = new TCanvas("c3","Exposition and Count on Equator",1200,3000); |
| 226 |
TCanvas *c4 = new TCanvas("c4","Pitch-Angle exposition and count",1200,3400); |
| 227 |
TCanvas *c5 = new TCanvas("c5","PitchAndLieveTimevsLatitude",1200,800); |
| 228 |
TCanvas *c6 = new TCanvas("c6","Diferense",1200,800); |
| 229 |
//TCanvas *c6 = new TCanvas("c6","CountEquator",1200,3000); |
| 230 |
|
| 231 |
//} |
| 232 |
|
| 233 |
/***************************************************************************************/ |
| 234 |
//Histogramms for diferences in results between calculating using DoTrack2 |
| 235 |
//and using axv[0], ayv[0] |
| 236 |
/**************************************************************************************/ |
| 237 |
|
| 238 |
//if(DiffHist){ |
| 239 |
|
| 240 |
TH1F *DifDoTr2ANdaxvAzim = new TH1F("DifDoTr2AndaxvAzim","Diference between Azimutal angles",360,-180,180); |
| 241 |
TH1F *DifDoTr2ANdaxvZenit = new TH1F("DifDoTr2AndaxvZenit","Diference between Zenit angles",360,-180,180); |
| 242 |
|
| 243 |
//} |
| 244 |
|
| 245 |
/*************************************************************************************/ |
| 246 |
//Histogramm for Pamela's angular efficiency |
| 247 |
/*************************************************************************************/ |
| 248 |
|
| 249 |
TH2F *PamAngEffhigh = new TH2F("PamAngEffhigh","Zenit Angle, deg; Azimutal Angle, deg",360,0,360,100,80,180); |
| 250 |
TH2F *PamAngEfflow = new TH2F("PamAngEfflow","Zenit Angle, deg; Azimutal Angle, deg",360,0,360,100,80,180); |
| 251 |
TCanvas *c7 = new TCanvas("c7","Pamela's angular efficiency", 1600,2000); |
| 252 |
|
| 253 |
/*************************************************************************************/ |
| 254 |
//Hystogramm to count how much time Pamela have each PitchAngle |
| 255 |
/*************************************************************************************/ |
| 256 |
|
| 257 |
TH1F *PitchTime = new TH1F("PitchTime","Time for Pitch",360,0,360); |
| 258 |
TCanvas *c8 = new TCanvas("c8","PitchTime",1200,600); |
| 259 |
|
| 260 |
TH1F *PitchExposure = new TH1F("PitchExposure","Time for Pitch",360,0,360); |
| 261 |
TCanvas *c9 = new TCanvas("c9","PitchTime",1200,600); |
| 262 |
|
| 263 |
/***************************************************************************************/ |
| 264 |
//Initialization histogramms for count of protons depends from Energy. Volodia. |
| 265 |
/***************************************************************************************/ |
| 266 |
/* |
| 267 |
TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.); |
| 268 |
proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV"); |
| 269 |
TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.); |
| 270 |
TH1F* binw= new TH1F("w","",80,-1.,3.); |
| 271 |
proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV"); |
| 272 |
proton1log->SetLineColor(kRed); |
| 273 |
TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800); |
| 274 |
// this is bins wide to calculate Flux |
| 275 |
Float_t x; |
| 276 |
Float_t xmin, xmax; |
| 277 |
Float_t binwide[80]; |
| 278 |
for (Int_t l=0 ; l<80; l++) { |
| 279 |
binw->Fill(-1.+0.05*l,pow(10.,0.05*(Float_t)(l+1)-1.)- pow(10.,0.05*(Float_t)(l)-1.)); |
| 280 |
} |
| 281 |
|
| 282 |
for (Int_t l=0 ; l<80; l++) { |
| 283 |
binwide[l]=binw->GetBinContent(l+1); |
| 284 |
} |
| 285 |
*/ |
| 286 |
/****************************************************************************/ |
| 287 |
//Geometrical Factor. Volodia |
| 288 |
/****************************************************************************/ |
| 289 |
/* |
| 290 |
TH1F* Geomfactor = new TH1F("G","G",1000,0.1,500.1); |
| 291 |
TH1F* Geomfactorlog= new TH1F("Glog","Glog",80,-1.,3.); |
| 292 |
TH1F* Geomfactorlogelectron= new TH1F("Gloge","Gloge",80,-1.,3.); |
| 293 |
for (Int_t l=0 ; l <80 ;l++) { |
| 294 |
x=pow(10.,(Float_t) 0.05*l-1.); |
| 295 |
Geomfactorlogelectron->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083))); |
| 296 |
x=ConvertT2R(x,0.938,1., 1.); |
| 297 |
Geomfactorlog->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083))); |
| 298 |
}//geomfactor for linear scale |
| 299 |
for (Int_t l=0 ; l <1000 ;l++) { |
| 300 |
x=(Float_t) 0.5*l+0.1; |
| 301 |
Geomfactor->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083))); |
| 302 |
} |
| 303 |
*/ |
| 304 |
/****************************************************************************/ |
| 305 |
//General loop |
| 306 |
/****************************************************************************/ |
| 307 |
|
| 308 |
for(ULong_t i=0; i<nevents; i++){ |
| 309 |
|
| 310 |
T->GetEntry(i); |
| 311 |
|
| 312 |
if(pam_events->GetTrkLevel2()->GetNTracks()==1){ |
| 313 |
|
| 314 |
/*************************************************************************/ |
| 315 |
//Getting general parameters |
| 316 |
/*************************************************************************/ |
| 317 |
|
| 318 |
track = pam_events->GetTrack(0); |
| 319 |
Int_t M0DE = pam_events->GetOrbitalInfo()->mode; //0 is zero here |
| 320 |
Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown; |
| 321 |
Float_t By = pam_events->GetOrbitalInfo()->Beast; |
| 322 |
Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth; |
| 323 |
Float_t Babs = pam_events->GetOrbitalInfo()->Babs; |
| 324 |
Float_t L = pam_events->GetOrbitalInfo()->L; |
| 325 |
Float_t rigev = 1/track->GetTrkTrack()->GetDeflection(); |
| 326 |
|
| 327 |
/*************************************************************************/ |
| 328 |
//Track selection |
| 329 |
/*************************************************************************/ |
| 330 |
|
| 331 |
ntr = 0; |
| 332 |
if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=500)) ntr=-1; |
| 333 |
for(Int_t ii=1;ii<=12;ii++){ |
| 334 |
if ((track->GetToFTrack()->beta[ii]<0) || (track->GetToFTrack()->beta[ii]>99)) ntr=-1; |
| 335 |
} |
| 336 |
if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1; |
| 337 |
if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1; |
| 338 |
if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1; |
| 339 |
//if(rigev<0.600) ntr=-1; |
| 340 |
|
| 341 |
if(ntr!=-1 || PamExp){ |
| 342 |
|
| 343 |
/*******************************************************************************/ |
| 344 |
//Transit reference system |
| 345 |
/*******************************************************************************/ |
| 346 |
|
| 347 |
TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(pam_events->GetOrbitalInfo()->q0,pam_events->GetOrbitalInfo()->q1,pam_events->GetOrbitalInfo()->q2,pam_events->GetOrbitalInfo()->q3),pam_events->GetOrbitalInfo()->absTime); |
| 348 |
TMatrixD Dij = PO->GreenwichtoGEO(pam_events->GetOrbitalInfo()->lat,pam_events->GetOrbitalInfo()->lon,Fij); |
| 349 |
TMatrixD Iij = PO->ColPermutation(Dij); |
| 350 |
|
| 351 |
/*****************************************************************************/ |
| 352 |
//Calculate Zenit and Azimutal angle and vector of particle flight in Pamela |
| 353 |
//using axv[0] and ayv[0] variables |
| 354 |
/*****************************************************************************/ |
| 355 |
|
| 356 |
if(!DoTr || DiffHist){ |
| 357 |
|
| 358 |
Double_t Aaxv = TMath::Abs(track->GetTrkTrack()->axv[0])*TMath::DegToRad(); |
| 359 |
Double_t Aayv = TMath::Abs(track->GetTrkTrack()->ayv[0])*TMath::DegToRad(); |
| 360 |
PamZenit = TMath::RadToDeg()*asin(sqrt(pow(sin(Aayv), 2) + pow(sin(Aaxv), 2))); |
| 361 |
|
| 362 |
Double_t axv = -track->GetTrkTrack()->axv[0] * TMath::DegToRad(); |
| 363 |
Double_t ayv = -track->GetTrkTrack()->ayv[0] * TMath::DegToRad(); |
| 364 |
Double_t angle = atan(sin(TMath::Abs(ayv))/sin(TMath::Abs(axv))) * TMath::RadToDeg(); |
| 365 |
|
| 366 |
PamAzim = 360. - angle; |
| 367 |
if(axv>=0 && ayv >=0) PamAzim = angle; |
| 368 |
if(axv<0 && ayv >0) PamAzim = 180. - angle; |
| 369 |
if(axv<0 && ayv <0) PamAzim = 180. + angle; |
| 370 |
|
| 371 |
PamAzim = PamAzim * TMath::DegToRad(); |
| 372 |
PamZenit = (180 - PamZenit) * TMath::DegToRad(); |
| 373 |
Px = sin(axv); |
| 374 |
Py = sin(ayv); |
| 375 |
Pz = cos(PamZenit); |
| 376 |
|
| 377 |
|
| 378 |
} |
| 379 |
|
| 380 |
/*****************************************************************************/ |
| 381 |
//Calculate Zenit and Azimutal angle and vector of particle flight in Pamela |
| 382 |
//using DoTrack2 procedure |
| 383 |
/*****************************************************************************/ |
| 384 |
|
| 385 |
if(DoTr){ |
| 386 |
|
| 387 |
track->GetTrkTrack()->DoTrack2(tr); |
| 388 |
Double_t E11x = tr->x[0]; |
| 389 |
Double_t E11y = tr->y[0]; |
| 390 |
Double_t E11z = zin[0]; |
| 391 |
Double_t E22x = tr->x[3]; |
| 392 |
Double_t E22y = tr->y[3]; |
| 393 |
Double_t E22z = zin[3]; |
| 394 |
Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2)); |
| 395 |
MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x)); |
| 396 |
MyAzim = MyAzim;//-180; |
| 397 |
MyZenit = TMath::RadToDeg()*acos((E22z-E11z)/norm); |
| 398 |
if(E22x-E11x>=0 && E22y-E11y <0) MyAzim = 360. - MyAzim; |
| 399 |
if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim; |
| 400 |
if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim; |
| 401 |
if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim; |
| 402 |
Px = (E22x-E11x)/norm; |
| 403 |
Py = (E22y-E11y)/norm; |
| 404 |
Pz = (E22z-E11z)/norm; |
| 405 |
|
| 406 |
} |
| 407 |
|
| 408 |
/***************************************************************************************/ |
| 409 |
//Fill histogramms for diferences in results between calculating using DoTrack2 |
| 410 |
//and using axv[0], ayv[0] |
| 411 |
/**************************************************************************************/ |
| 412 |
|
| 413 |
if(DiffHist){ |
| 414 |
|
| 415 |
DifDoTr2ANdaxvZenit->Fill(MyZenit-(TMath::RadToDeg()*PamZenit)); |
| 416 |
DifDoTr2ANdaxvAzim->Fill(MyAzim-(TMath::RadToDeg()*PamAzim)); |
| 417 |
|
| 418 |
} |
| 419 |
|
| 420 |
if(PamEff||AnglHist||PamExp){ |
| 421 |
|
| 422 |
/***********************************************************************************/ |
| 423 |
//Transit vector in Pamela reference system to GEO |
| 424 |
/***********************************************************************************/ |
| 425 |
|
| 426 |
TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz); |
| 427 |
Double_t A1 = Iij(0,2); |
| 428 |
Double_t A2 = Iij(1,2); |
| 429 |
Double_t A3 = Iij(2,2); |
| 430 |
|
| 431 |
/*******************************************************************************/ |
| 432 |
//Calculating angles (Pitch, Pamela's main axis etc.) |
| 433 |
/*******************************************************************************/ |
| 434 |
|
| 435 |
Double_t SB = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz); // Pitch angle |
| 436 |
Double_t SA = PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz); // Angle between Pamela's main axiz and B |
| 437 |
Double_t SC = PO->GetPitchAngle(1,0,0,Bx,By,Bz); // Angle between direction to xenit and B |
| 438 |
Double_t ZC = PO->GetPitchAngle(Px,Py,Pz,0,0,1); // Zenit Angle of Particle flight |
| 439 |
|
| 440 |
/******************************************************************************/ |
| 441 |
//Proton selection. Volodia. |
| 442 |
/******************************************************************************/ |
| 443 |
/* |
| 444 |
if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) { |
| 445 |
Int_t detr = 0; |
| 446 |
//if((rigev < -0.8 && betaev >0.8)||(rigev>-0.8 && rigev < 0.&& betaev > 0.8) ) //electron2->Fill(log10(fabs(rigev))); |
| 447 |
if(detr < (15*exp(-1.6*(fabs(rigev)-0.5))+4) && rigev >0.5 || (detr < 30 && fabs(rigev) <0.5)) { |
| 448 |
Float_t Ekin=ConvertR2T(rigev,0.938,1., 1.); |
| 449 |
proton1log->Fill(log10(Ekin)); |
| 450 |
} |
| 451 |
} |
| 452 |
*/ |
| 453 |
/************/ |
| 454 |
//Exposure |
| 455 |
/************/ |
| 456 |
|
| 457 |
//for (Int_t itime=0; itime<16;itime++){ |
| 458 |
// if(14.9/L/L > itime && 14.9/L/L < (itime+1) & Babs >0.24) exposure[itime]=exposure[itime]+livetime/1000.;} |
| 459 |
|
| 460 |
/***********************************************************************/ |
| 461 |
//Polar area. To calculate angular efficiencÐÐy of registration |
| 462 |
/***********************************************************************/ |
| 463 |
|
| 464 |
if(PamExp){ |
| 465 |
|
| 466 |
Float_t AnglTime = 0.16*pam_events->GetTrigLevel2()->dltime[0]+0.01*pam_events->GetTrigLevel2()->dltime[1]; |
| 467 |
for(Int_t i=0;i<=360;i++){ |
| 468 |
for(Int_t j=0;j<=100;j++){ |
| 469 |
if(hhigh->GetBin(i,j)>0){ |
| 470 |
Double_t pamZenit=(j+80-0.5)*TMath::DegToRad(); |
| 471 |
Double_t pamAzim=(i*-0.5)*TMath::DegToRad(); |
| 472 |
Double_t _z = cos(pamZenit); |
| 473 |
Double_t _x = sin(pamZenit)*cos(pamAzim); |
| 474 |
Double_t _y = sin(pamZenit)*sin(pamAzim); |
| 475 |
TMatrixD Uij = PO->PamelatoGEO(Iij,_x,_y,_z); |
| 476 |
Double_t SP = PO->GetPitchAngle(Uij(0,0),Uij(1,0),Uij(2,0),Bx,By,Bz); |
| 477 |
PitchExposure->Fill(SP,AnglTime); |
| 478 |
} |
| 479 |
} |
| 480 |
} |
| 481 |
|
| 482 |
} |
| 483 |
|
| 484 |
if(PamEff){ |
| 485 |
|
| 486 |
if(L>6){ |
| 487 |
Float_t Mass = 0; |
| 488 |
Int_t detr = 0; |
| 489 |
Float_t Ekin = 0; |
| 490 |
Float_t betaev = 0.25*(track->GetToFTrack()->beta[0]+track->GetToFTrack()->beta[1]+track->GetToFTrack()->beta[2]+track->GetToFTrack()->beta[3]); |
| 491 |
if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) { |
| 492 |
if((rigev < -0.8 && betaev >0.8)||(rigev>-0.8 && rigev < 0.&& betaev > 0.8) ) { |
| 493 |
Ekin=ConvertR2T(rigev,0.511,-1.); |
| 494 |
if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit); |
| 495 |
if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit); |
| 496 |
}//electron2->Fill(log10(fabs(rigev))); |
| 497 |
|
| 498 |
if(detr < (15*exp(-1.6*(fabs(rigev)-0.5))+4) && rigev >0.5 || (detr < 30 && fabs(rigev) <0.5)) { |
| 499 |
Ekin=ConvertR2T(rigev,0.938,-1.); |
| 500 |
if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit); |
| 501 |
if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit); |
| 502 |
} |
| 503 |
} |
| 504 |
} |
| 505 |
|
| 506 |
} |
| 507 |
|
| 508 |
/***********************************************************************/ |
| 509 |
//Filling hystogram to count how much time Pamela have each Pitch-Angle |
| 510 |
/***********************************************************************/ |
| 511 |
|
| 512 |
if(PamAngTime) { |
| 513 |
|
| 514 |
Float_t AnglTime = 0.16*pam_events->GetTrigLevel2()->dltime[0]+0.01*pam_events->GetTrigLevel2()->dltime[1]; |
| 515 |
PitchTime->Fill(SA,AnglTime); |
| 516 |
|
| 517 |
} |
| 518 |
|
| 519 |
if(AnglHist){ |
| 520 |
|
| 521 |
/***********************************************************************/ |
| 522 |
//Filling angular histogramms for Brazil |
| 523 |
/***********************************************************************/ |
| 524 |
|
| 525 |
if((Babs<0.19)&&(L<1.75)){ |
| 526 |
Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0]; |
| 527 |
PitchExpositionBrazil->Fill(SB,Argv); |
| 528 |
MainAxisPamelaExpositionBrazil->Fill(SA,Argv); |
| 529 |
PitchExpositionNoOrientationBrazil->Fill(SC,Argv); |
| 530 |
PitchCountBrazil->Fill(SB); |
| 531 |
MainAxisPamelaCountBrazil->Fill(SA); |
| 532 |
PitchCountNoOrientationBrazil->Fill(SC); |
| 533 |
ZenitAngle->Fill(ZC); |
| 534 |
PitchvsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,SA); |
| 535 |
LieveTimevsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,Argv); |
| 536 |
} |
| 537 |
|
| 538 |
/***********************************************************************/ |
| 539 |
//Filling angulars histogramms for areas without Brazil |
| 540 |
/***********************************************************************/ |
| 541 |
|
| 542 |
if(Babs>0.24){ |
| 543 |
Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0]; |
| 544 |
PitchExpositionWithoutBrazil->Fill(SB,Argv); |
| 545 |
MainAxisPamelaExpositionWithoutBrazil->Fill(SA,Argv); |
| 546 |
PitchExpositionNoOrientationWithoutBrazil->Fill(SC,Argv); |
| 547 |
PitchCountWithoutBrazil->Fill(SB); |
| 548 |
MainAxisPamelaCountWithoutBrazil->Fill(SA); |
| 549 |
PitchCountNoOrientationWithoutBrazil->Fill(SC); |
| 550 |
ZenitAngle->Fill(ZC); |
| 551 |
|
| 552 |
/***************************************************************/ |
| 553 |
//Filling angulars histogramm for equator |
| 554 |
/***************************************************************/ |
| 555 |
|
| 556 |
if(L<1.2){ |
| 557 |
Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0]; |
| 558 |
PitchExpositionEquator->Fill(SB,Argv); |
| 559 |
MainAxisPamelaExpositionEquator->Fill(SA,Argv); |
| 560 |
PitchExpositionNoOrientationEquator->Fill(SC,Argv); |
| 561 |
PitchCountEquator->Fill(SB); |
| 562 |
MainAxisPamelaCountEquator->Fill(SA); |
| 563 |
PitchCountNoOrientationEquator->Fill(SC); |
| 564 |
ZenitAngle->Fill(ZC); |
| 565 |
} |
| 566 |
} |
| 567 |
|
| 568 |
} // if(AnglHist) |
| 569 |
|
| 570 |
} // if(PamEff||AnglHist) |
| 571 |
|
| 572 |
} // if ntr!=-1; |
| 573 |
|
| 574 |
} // if GetNTrack==1; |
| 575 |
|
| 576 |
} //general loop |
| 577 |
|
| 578 |
Double_t BinContent; |
| 579 |
|
| 580 |
/*********************************************************/ |
| 581 |
//Divide canvases for angular histogramms |
| 582 |
/*********************************************************/ |
| 583 |
|
| 584 |
if(AnglHist){ |
| 585 |
|
| 586 |
c1->Divide(1,6); |
| 587 |
c2->Divide(1,6); |
| 588 |
c3->Divide(1,6); |
| 589 |
c4->Divide(1,7); |
| 590 |
c5->Divide(1,2); |
| 591 |
|
| 592 |
} |
| 593 |
|
| 594 |
/***************************************************************************************/ |
| 595 |
//Canvas for histogramms for diferences in results between calculating using DoTrack2 |
| 596 |
//and using axv[0], ayv[0] |
| 597 |
/**************************************************************************************/ |
| 598 |
|
| 599 |
if(DiffHist){ |
| 600 |
|
| 601 |
c6->Divide(1,2); |
| 602 |
|
| 603 |
} |
| 604 |
|
| 605 |
/**************************************************************************************/ |
| 606 |
//Canvas for Pamela's efficiency hystogramm |
| 607 |
/**************************************************************************************/ |
| 608 |
|
| 609 |
if(PamEff){ |
| 610 |
|
| 611 |
c7->Divide(1,2); |
| 612 |
|
| 613 |
} |
| 614 |
|
| 615 |
/**************************************************************************************/ |
| 616 |
//Draw Angular Histogramms |
| 617 |
/**************************************************************************************/ |
| 618 |
|
| 619 |
if(AnglHist){ |
| 620 |
|
| 621 |
c1->cd(1); |
| 622 |
MainAxisPamelaExpositionBrazil->Draw(); |
| 623 |
c1->cd(2); |
| 624 |
MainAxisPamelaCountBrazil->Draw(); |
| 625 |
c1->cd(3); |
| 626 |
PitchExpositionNoOrientationBrazil->Draw(); |
| 627 |
c1->cd(4); |
| 628 |
PitchCountNoOrientationBrazil->Draw(); |
| 629 |
c1->cd(5); |
| 630 |
DiferenceBetweenNoOrientationAndWithOrientationBrazil->Divide(PitchExpositionNoOrientationBrazil,MainAxisPamelaExpositionBrazil); |
| 631 |
DiferenceBetweenNoOrientationAndWithOrientationBrazil->Draw(); |
| 632 |
c1->cd(6); |
| 633 |
for(Int_t iu=0;iu<360;iu++){ |
| 634 |
BinContent = DiferenceBetweenNoOrientationAndWithOrientationBrazil->GetBinContent(iu); |
| 635 |
DifferenceCountInBrazil->Fill(BinContent); |
| 636 |
} |
| 637 |
DifferenceCountInBrazil->Draw(); |
| 638 |
//DifferenceCountInBrazil->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.root"); |
| 639 |
c1->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.jpg"); |
| 640 |
|
| 641 |
c2->cd(1); |
| 642 |
MainAxisPamelaExpositionWithoutBrazil->Draw(); |
| 643 |
c2->cd(2); |
| 644 |
MainAxisPamelaCountWithoutBrazil->Draw(); |
| 645 |
c2->cd(3); |
| 646 |
PitchExpositionNoOrientationWithoutBrazil->Draw(); |
| 647 |
c2->cd(4); |
| 648 |
PitchCountNoOrientationWithoutBrazil->Draw(); |
| 649 |
c2->cd(5); |
| 650 |
DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Divide(PitchExpositionNoOrientationWithoutBrazil,MainAxisPamelaExpositionWithoutBrazil); |
| 651 |
DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Draw(); |
| 652 |
c2->cd(6); |
| 653 |
for(Int_t iu=0;iu<360;iu++){ |
| 654 |
BinContent = DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->GetBinContent(iu); |
| 655 |
DifferenceCountWithoutBrazil->Fill(BinContent); |
| 656 |
} |
| 657 |
DifferenceCountWithoutBrazil->Draw(); |
| 658 |
c2->SaveAs((TString)outdir+"PIC003A_WithoutBrazilExpositionAndCount.jpg"); |
| 659 |
|
| 660 |
c3->cd(1); |
| 661 |
MainAxisPamelaExpositionEquator->Draw(); |
| 662 |
c3->cd(2); |
| 663 |
MainAxisPamelaCountEquator->Draw(); |
| 664 |
c3->cd(3); |
| 665 |
PitchExpositionNoOrientationEquator->Draw(); |
| 666 |
c3->cd(4); |
| 667 |
PitchCountNoOrientationEquator->Draw(); |
| 668 |
c3->cd(5); |
| 669 |
DiferenceBetweenNoOrientationAndWithOrientationEquator->Divide(PitchExpositionNoOrientationEquator,MainAxisPamelaExpositionEquator); |
| 670 |
DiferenceBetweenNoOrientationAndWithOrientationEquator->Draw(); |
| 671 |
c3->cd(6); |
| 672 |
for(Int_t iu=0;iu<360;iu++){ |
| 673 |
BinContent = DiferenceBetweenNoOrientationAndWithOrientationEquator->GetBinContent(iu); |
| 674 |
DifferenceCountEquator->Fill(BinContent); |
| 675 |
} |
| 676 |
DifferenceCountEquator->Draw(); |
| 677 |
c3->SaveAs((TString)outdir+"PIC005A_EquatorExpositionAndCount.jpg"); |
| 678 |
|
| 679 |
c4->cd(1); |
| 680 |
PitchExpositionBrazil->Draw(); |
| 681 |
c4->cd(2); |
| 682 |
PitchCountBrazil->Draw(); |
| 683 |
c4->cd(3); |
| 684 |
PitchExpositionWithoutBrazil->Draw(); |
| 685 |
c4->cd(4); |
| 686 |
PitchCountWithoutBrazil->Draw(); |
| 687 |
c4->cd(5); |
| 688 |
PitchExpositionEquator->Draw(); |
| 689 |
c4->cd(6); |
| 690 |
PitchCountEquator->Draw(); |
| 691 |
c4->cd(7); |
| 692 |
ZenitAngle->Draw(); |
| 693 |
c4->SaveAs((TString)outdir+"PIC006A_PitchCountsAndExpositions.jpg"); |
| 694 |
|
| 695 |
c5->cd(1); |
| 696 |
PitchvsLatInBrazil->Draw(); |
| 697 |
c5->cd(2); |
| 698 |
LieveTimevsLatInBrazil->Draw(); |
| 699 |
c5->SaveAs((TString)outdir+"PIC006A_PitchvsLotitudeBrazil.jpg"); |
| 700 |
|
| 701 |
} |
| 702 |
|
| 703 |
/***************************************************************************************/ |
| 704 |
//Draw hystogramms for diferences in results between calculating using DoTrack2 |
| 705 |
//and using axv[0], ayv[0] |
| 706 |
/**************************************************************************************/ |
| 707 |
|
| 708 |
if(DiffHist){ |
| 709 |
|
| 710 |
c6->cd(1); |
| 711 |
DifDoTr2ANdaxvAzim->Draw(); |
| 712 |
c6->cd(2); |
| 713 |
DifDoTr2ANdaxvZenit->Draw(); |
| 714 |
c6->SaveAs((TString)outdir+"PIC001Diference.jpg"); |
| 715 |
|
| 716 |
} |
| 717 |
|
| 718 |
/***********************************************************************/ |
| 719 |
//Draw hystogram to count how much time Pamela have each Pitch-Angle |
| 720 |
/***********************************************************************/ |
| 721 |
|
| 722 |
if(PamAngTime) { |
| 723 |
|
| 724 |
c8->cd(); |
| 725 |
PitchTime->Draw(); |
| 726 |
c8->SaveAs((TString)outdir+"PIC002PamAngTime.jpg"); |
| 727 |
|
| 728 |
} |
| 729 |
|
| 730 |
if(PamExp){ |
| 731 |
|
| 732 |
c9->cd(); |
| 733 |
PitchExposure->Draw(); |
| 734 |
c9->SaveAs((TString)outdir+"PIC003PamAngTime.jpg"); |
| 735 |
|
| 736 |
} |
| 737 |
|
| 738 |
/**************************************************************************************/ |
| 739 |
//Draw Hystogramm for Pamela's angular efficiency |
| 740 |
/**************************************************************************************/ |
| 741 |
gStyle->SetPalette(1); |
| 742 |
|
| 743 |
if(PamEff){ |
| 744 |
c7->cd(1); |
| 745 |
PamAngEffhigh->Draw("colz"); |
| 746 |
TFile fh((TString)outdir+"PamAngEfficiencyHighEnergy.root","recreate"); |
| 747 |
PamAngEffhigh->Write(); |
| 748 |
fh.Close(); |
| 749 |
//PamAngEffhigh->SaveAs((TString)outdir+"PamAngEfficiencyHighEnergy.root"); |
| 750 |
c7->cd(2); |
| 751 |
PamAngEfflow->Draw("colz"); |
| 752 |
TFile fl((TString)outdir+"PamAngEfficiencyLowEnergy.root","recreate"); |
| 753 |
PamAngEffhigh->Write(); |
| 754 |
fl.Close(); |
| 755 |
//PamAngEfflow->SaveAs((TString)outdir+"PamAngEfficiencyLowEnergy.root"); |
| 756 |
c7->SaveAs((TString)outdir+"PIC001PamAngEfficiencyHighEnergy.jpg"); |
| 757 |
|
| 758 |
} |
| 759 |
|
| 760 |
/**************************************************************************************/ |
| 761 |
//Draw hystogramms for protons. Volodia |
| 762 |
/**************************************************************************************/ |
| 763 |
|
| 764 |
/* |
| 765 |
plogCanvas->cd(); |
| 766 |
//plogCanvas->SetLogx(); |
| 767 |
plogCanvas->SetLogy(); |
| 768 |
plogCanvas->SetGrid(); |
| 769 |
//BinLogX(proton1); |
| 770 |
Float_t Energybin[80],Fluxbin[80],dE[80],dFlux[80]; |
| 771 |
|
| 772 |
for (Int_t i=0;i<80;i++){ |
| 773 |
dFlux[i]=proton1log->GetBinContent(i+1); |
| 774 |
dFlux[i]=1./sqrt(dFlux[i]); |
| 775 |
// bintime3eq[i]=time3eq->GetBinContent(i+1); |
| 776 |
} |
| 777 |
|
| 778 |
proton1log->Divide(proton1log,binw); |
| 779 |
proton1log->Divide(proton1log,Geomfactorlog); |
| 780 |
proton1log->Scale(0.001); //MeV ->GeV |
| 781 |
proton1log->Scale(10000.); //cm2 -> m2 |
| 782 |
//proton1log->Scale(1./exposure[0]); |
| 783 |
|
| 784 |
for (Int_t i=0;i<80;i++){ |
| 785 |
Fluxbin[i]=proton1log->GetBinContent(i+1); |
| 786 |
dFlux[i]=Fluxbin[i]*dFlux[i]; |
| 787 |
Energybin[i]=pow(10.,(Float_t)i*0.05+0.025-1); |
| 788 |
Energybin[i+1]=pow(10.,(Float_t)(i+1)*0.05+0.025-1.); |
| 789 |
dE[i]=(Energybin[i+1]-Energybin[i])/2.; |
| 790 |
// bintime3eq[i]=time3eq->GetBinContent(i+1); |
| 791 |
//cout<<i<<" "<<bintime3[i]<<" "<<bintime3eq[i]<<endl; |
| 792 |
//flux_out<<Energybin[i]<<" "<<Fluxbin[i]<<" "<<dE[i]<<" "<<dFlux[i]<<endl; |
| 793 |
} |
| 794 |
|
| 795 |
proton1log->Draw(""); |
| 796 |
|
| 797 |
*/ |
| 798 |
} |
Parent Directory
| 
Revision Log