| 1 |
pam-mep |
1.1 |
/* |
| 2 |
|
|
* PitchAngleExpositionAction.cpp |
| 3 |
|
|
* |
| 4 |
|
|
* Created on: 04-april-2010 |
| 5 |
|
|
* Author: Vitaly Malakhov |
| 6 |
|
|
*/ |
| 7 |
|
|
|
| 8 |
|
|
/*! @file PitchAngleExpositionAction.cpp The PitchAngleExpositionAction class implementation file. */ |
| 9 |
|
|
|
| 10 |
|
|
#include "PitchAngleExpositionAction.h" |
| 11 |
|
|
#include "TMatrixD.h" |
| 12 |
|
|
#include "OrientationInfo.h" |
| 13 |
|
|
|
| 14 |
|
|
#define XBINS_F 3600 |
| 15 |
|
|
#define YBINS_F 400 |
| 16 |
|
|
|
| 17 |
|
|
PitchAngleExpositionAction::PitchAngleExpositionAction(const char *actionName, TString outFileBase, TString mode, |
| 18 |
|
|
bool outRoot, bool outText, TString title, bool CalcPitchExpo, UInt_t nPitchBins, Float_t lowerPitch, Float_t upperPitch) : |
| 19 |
|
|
Histo2DAction<Float_t> (actionName, title, outFileBase, mode, outRoot, outText), _calcPitchExpo(CalcPitchExpo), _nPitchBins(nPitchBins), _lowerPitch(lowerPitch), _upperPitch(upperPitch) { |
| 20 |
|
|
|
| 21 |
|
|
} |
| 22 |
|
|
|
| 23 |
|
|
void PitchAngleExpositionAction::OnGood(PamLevel2 *event) { |
| 24 |
|
|
Float_t Bx = -event->GetOrbitalInfo()->Bdown; |
| 25 |
|
|
Float_t By = event->GetOrbitalInfo()->Beast; |
| 26 |
|
|
Float_t Bz = event->GetOrbitalInfo()->Bnorth; |
| 27 |
|
|
OrientationInfo* PO = new OrientationInfo(); |
| 28 |
|
|
TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(event->GetOrbitalInfo()->q0,event->GetOrbitalInfo()->q1,event->GetOrbitalInfo()->q2,event->GetOrbitalInfo()->q3),event->GetOrbitalInfo()->absTime); |
| 29 |
|
|
TMatrixD Dij = PO->GreenwichtoGEO(event->GetOrbitalInfo()->lat,event->GetOrbitalInfo()->lon,Fij); |
| 30 |
|
|
TMatrixD Iij = PO->ColPermutation(Dij); |
| 31 |
|
|
TMatrixD Sij = Iij; |
| 32 |
|
|
Sij.Invert(); |
| 33 |
|
|
TMatrixD Nij = PO->PamelatoGEO(Sij,Bx,By,Bz); |
| 34 |
|
|
Float_t BZen = acos(Nij(2,0)) * TMath::RadToDeg(); |
| 35 |
|
|
Float_t angle = 0; |
| 36 |
|
|
if(Nij(0,0)!=0)angle = atan(Nij(1,0)/Nij(0,0)) * TMath::RadToDeg();else angle = 90.; |
| 37 |
|
|
Float_t BAz = 360. + angle; |
| 38 |
|
|
if(Nij(0,0) >=0 && Nij(1,0) >=0) BAz = angle; |
| 39 |
|
|
if(Nij(0,0)<0) BAz = 180. + angle; |
| 40 |
|
|
Float_t weight = event->GetTrigLevel2()->dltime[0] * 0.16 /1000.; |
| 41 |
|
|
Fill(BAz, BZen, weight); |
| 42 |
|
|
PO->Delete(); |
| 43 |
|
|
} |
| 44 |
|
|
|
| 45 |
|
|
void PitchAngleExpositionAction::Finalize() { |
| 46 |
|
|
|
| 47 |
|
|
Histo2DAction<Float_t>::Finalize(); |
| 48 |
|
|
|
| 49 |
|
|
if(_calcPitchExpo){ |
| 50 |
|
|
TH2F tmphisto; |
| 51 |
|
|
tmphisto.SetBins(XBINS_F,0,360,YBINS_F,140,180); |
| 52 |
|
|
TString tmpname = GetRootHisto()->GetName(); |
| 53 |
|
|
Float_t delta = (_upperPitch - _lowerPitch)/_nPitchBins; |
| 54 |
|
|
for(UInt_t i=0; i<_nPitchBins; i++){ |
| 55 |
|
|
_pitchAngleHisto.push_back(tmphisto); |
| 56 |
|
|
TString tmp = tmpname; |
| 57 |
|
|
TString tmp2 = ""; |
| 58 |
|
|
tmp+=i*delta; |
| 59 |
|
|
tmp+="_"; |
| 60 |
|
|
tmp+=(i+1)*delta; |
| 61 |
|
|
tmp+="deg"; |
| 62 |
|
|
for(Int_t j = 0; j<tmp.Sizeof()-1; j++){ |
| 63 |
|
|
if(tmp[j] == '.') tmp[j] = 'p'; |
| 64 |
|
|
if(tmp[j] == ' '){ |
| 65 |
|
|
tmp.Replace(j,1,NULL); |
| 66 |
|
|
j--; |
| 67 |
|
|
} |
| 68 |
|
|
} |
| 69 |
|
|
cout<<"size is "<<_pitchAngleHisto.size()<<"\t"<<tmp<<endl; |
| 70 |
|
|
_pitchAngleHisto[i].SetName(tmp); |
| 71 |
|
|
_pitchAngleHisto[i].SetTitle(tmp); |
| 72 |
|
|
cout<<"PitchHisto["<<i<<"] name is "<<_pitchAngleHisto[i].GetName()<<endl; |
| 73 |
|
|
} |
| 74 |
|
|
|
| 75 |
|
|
Float_t BinWidthX = GetRootHisto()->GetXaxis()->GetXmax()/GetRootHisto()->GetNbinsX(); |
| 76 |
|
|
Float_t BinWidthY = GetRootHisto()->GetYaxis()->GetXmax()/GetRootHisto()->GetNbinsY(); |
| 77 |
|
|
Float_t tetaB = 0.; //Zenith angle of Magnetic field vector B in PAMELA reference frame |
| 78 |
|
|
Float_t fiB = 0.; //Azimuth angle of Magnetic field vector B in PAMELA reference frame |
| 79 |
|
|
Float_t Psi = 0.; //Pitch-Angle; |
| 80 |
|
|
Float_t Ksi = 0.; //Azimuth-angle according pitch-angle Psi |
| 81 |
|
|
|
| 82 |
|
|
Bool_t pth[XBINS_F][YBINS_F]; |
| 83 |
|
|
for(Int_t i=1;i<=GetRootHisto()->GetNbinsX();i++){ |
| 84 |
|
|
cout<<i<<" of "<<GetRootHisto()->GetNbinsX()<<endl; |
| 85 |
|
|
for(Int_t j=1;j<=GetRootHisto()->GetNbinsY();j++){ |
| 86 |
|
|
if(GetRootHisto()->GetBinContent(i,j)!=0){ |
| 87 |
|
|
cout<<"j = "<<j<<endl; |
| 88 |
|
|
fiB = TMath::DegToRad()*(i*BinWidthX+BinWidthX/2); |
| 89 |
|
|
tetaB = TMath::DegToRad()*(j*BinWidthY+BinWidthY/2); |
| 90 |
|
|
for(UInt_t iPsi=0;iPsi<_nPitchBins;iPsi++){ //We suggest here that Pitch-angel bin size is 5 degrees; |
| 91 |
|
|
for(UInt_t ip=0;ip<XBINS_F;ip++){ |
| 92 |
|
|
for(UInt_t jp=0;jp<YBINS_F;jp++) pth[ip][jp]=false; |
| 93 |
|
|
} |
| 94 |
|
|
cout<<"iPsi = "<<iPsi<<endl; |
| 95 |
|
|
for(UInt_t jPsi=0;jPsi<delta*10;jPsi++){ //cycle inside one pitch-angle bin; |
| 96 |
|
|
for(UInt_t iKsi=0;iKsi<XBINS_F;iKsi++){ //cycle over azimuth angle according this Pitch-Angle; |
| 97 |
|
|
Psi=TMath::DegToRad()*(iPsi*5+jPsi*0.1); |
| 98 |
|
|
Ksi=TMath::DegToRad()*(iKsi*0.1); |
| 99 |
|
|
Float_t z_ = 0; |
| 100 |
|
|
Float_t x_ = sin(Psi)*cos(Ksi); |
| 101 |
|
|
Float_t y_ = sin(Psi)*sin(Ksi); |
| 102 |
|
|
Float_t x1 = x_*cos(tetaB)+z_*sin(tetaB); |
| 103 |
|
|
Float_t z1 = -x_*sin(tetaB)+z_*cos(tetaB); |
| 104 |
|
|
Float_t y1 = y_; |
| 105 |
|
|
Float_t x2 = x1*cos(fiB)+y1*sin(fiB); |
| 106 |
|
|
Float_t y2 = -x1*sin(fiB)+y1*cos(fiB); |
| 107 |
|
|
Float_t z2 = z1; |
| 108 |
|
|
Float_t a2 = pow(sin(Psi),2); |
| 109 |
|
|
Float_t b2 = 2-2*cos(Psi); |
| 110 |
|
|
Float_t l = 1 - sqrt(b2-a2); |
| 111 |
|
|
Float_t xm = l*sin(tetaB)*cos(fiB); |
| 112 |
|
|
Float_t ym = l*sin(tetaB)*sin(fiB); |
| 113 |
|
|
Float_t zm = l*cos(tetaB); |
| 114 |
|
|
Float_t x3 = x2 - xm; |
| 115 |
|
|
Float_t y3 = y2 - ym; |
| 116 |
|
|
Float_t z3 = z2 - zm; |
| 117 |
|
|
Float_t tetaDir = 0.; //tetaDir and fiDir is a direction for Pitch and Azimuth angle given in Instrumental reference frame |
| 118 |
|
|
if(z3!=0)tetaDir = TMath::RadToDeg()*atan(sqrt(pow(x3,2)+pow(y3,2))/z3); |
| 119 |
|
|
if(z3<0)tetaDir=180+tetaDir; |
| 120 |
|
|
if(z3==0)tetaDir=90.; |
| 121 |
|
|
tetaDir=180-tetaDir; |
| 122 |
|
|
Float_t fiDir = 0; |
| 123 |
|
|
if(x3!=0)fiDir = TMath::RadToDeg()*atan(TMath::Abs(y3)/TMath::Abs(x3)); |
| 124 |
|
|
if(x3<0 && y3>0) fiDir = 180. - fiDir; |
| 125 |
|
|
if(x3<0 && y3<0) fiDir+=180.; |
| 126 |
|
|
if(x3>0 && y3<0) fiDir = 360. - fiDir; |
| 127 |
|
|
if(x3==0 && y3>0) fiDir = 90.; |
| 128 |
|
|
if(x3==0 && y3<0) fiDir = 270.; |
| 129 |
|
|
if(x3<0 && y3==0) fiDir = 180.; |
| 130 |
|
|
|
| 131 |
|
|
if(!pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]){ |
| 132 |
|
|
_pitchAngleHisto[iPsi].Fill(fiDir,tetaDir,GetRootHisto()->GetBinContent(i,j)); |
| 133 |
|
|
pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]=true; |
| 134 |
|
|
}/* |
| 135 |
|
|
if(fiDir<=35.){ |
| 136 |
|
|
if(tetaDir>=158.8-0.0035*pow(fiDir,2)){ |
| 137 |
|
|
if(!pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]){ |
| 138 |
|
|
_pitchAngleHisto[iPsi].Fill(fiDir,tetaDir,GetRootHisto()->GetBinContent(i,j)); |
| 139 |
|
|
pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]=true; |
| 140 |
|
|
} |
| 141 |
|
|
} |
| 142 |
|
|
} |
| 143 |
|
|
if(fiDir>35. && fiDir<=140.){ |
| 144 |
|
|
if(tetaDir>=163.2-0.0035*pow(fiDir-90.,2)){ |
| 145 |
|
|
if(!pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]){ |
| 146 |
|
|
_pitchAngleHisto[iPsi].Fill(fiDir,tetaDir,GetRootHisto()->GetBinContent(i,j)); |
| 147 |
|
|
pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]=true; |
| 148 |
|
|
} |
| 149 |
|
|
} |
| 150 |
|
|
} |
| 151 |
|
|
if(fiDir>140. && fiDir<=220.){ |
| 152 |
|
|
if(tetaDir>=161.5-0.0035*pow(fiDir-180.,2)){ |
| 153 |
|
|
if(!pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]){ |
| 154 |
|
|
_pitchAngleHisto[iPsi].Fill(fiDir,tetaDir,GetRootHisto()->GetBinContent(i,j)); |
| 155 |
|
|
pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]=true; |
| 156 |
|
|
} |
| 157 |
|
|
} |
| 158 |
|
|
} |
| 159 |
|
|
if(fiDir>220. && fiDir<=320.){ |
| 160 |
|
|
if(tetaDir>=163.2-0.0035*pow(fiDir-270.,2)){ |
| 161 |
|
|
if(!pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]){ |
| 162 |
|
|
_pitchAngleHisto[iPsi].Fill(fiDir,tetaDir,GetRootHisto()->GetBinContent(i,j)); |
| 163 |
|
|
pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]=true; |
| 164 |
|
|
} |
| 165 |
|
|
} |
| 166 |
|
|
} |
| 167 |
|
|
if(fiDir>320. && fiDir<=360.){ |
| 168 |
|
|
if(tetaDir>=158.8-0.0035*pow(fiDir-360.,2)){ |
| 169 |
|
|
if(!(pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1])){ |
| 170 |
|
|
_pitchAngleHisto[iPsi].Fill(fiDir,tetaDir,GetRootHisto()->GetBinContent(i,j)); |
| 171 |
|
|
pth[(Int_t)(fiDir/0.1)+1][(Int_t)((tetaDir-140.)/0.1)+1]=true; |
| 172 |
|
|
} |
| 173 |
|
|
} |
| 174 |
|
|
}*/ |
| 175 |
|
|
} // end of loop over iKsi |
| 176 |
|
|
} // end of loop over jPsi |
| 177 |
|
|
} // end of loop over iPsi |
| 178 |
|
|
} // if(h->GetBinContent(i,j)!=0) |
| 179 |
|
|
} // end of loop over j |
| 180 |
|
|
} // end of loop over i |
| 181 |
|
|
TFile outRootFile((Histo2DAction<Float_t>::_outFileBase + ".root"), "UPDATE"); |
| 182 |
|
|
outRootFile.cd(); |
| 183 |
|
|
for (UInt_t i=0; i < _nPitchBins ;i++){ |
| 184 |
|
|
_pitchAngleHisto[i].Write(); |
| 185 |
|
|
cout<<i<<"\t"<<_pitchAngleHisto[i].GetName()<<endl; |
| 186 |
|
|
} |
| 187 |
|
|
outRootFile.Close(); |
| 188 |
|
|
|
| 189 |
|
|
cout<<"outFile writing"<<endl; |
| 190 |
|
|
} |
| 191 |
|
|
} |
Parent Directory
| 
Revision Log