tpc/generator/generate_distortion_map.C

0001
0002
0003 #include "AnnularFieldSim.h"
0004 #include "TTree.h" //this prevents a lazy binding issue and/or is a magic spell.
0005 #include "TCanvas.h" //this prevents a lazy binding issue and/or is a magic spell.
0006
0007 // cppcheck-suppress unknownMacro
0008 R__LOAD_LIBRARY(libfieldsim.so)
0009
0010 char field_string[200];
0011 char lookup_string[200];
0012
0013 AnnularFieldSim *SetupDefaultSphenixTpc(bool twinMe=false, bool useSpacecharge=true, float xshift=0, float yshift=0, float zshift=0);
0014 AnnularFieldSim *SetupDigitalCurrentSphenixTpc(bool twinMe=false, bool useSpacecharge=true);
0015 void TestSpotDistortion(AnnularFieldSim *t);
0016 void SurveyFiles(TFileCollection* filelist);
0017
0018
0019 void generate_distortion_map(const char *inputname, const char* gainName, const char *outputname, const char *ibfName, const char *primName, bool hasSpacecharge=true, bool isAdc=false, int nSteps=500, bool scanSteps=false, float xshift=0, float yshift=0, float zshift=0){
0020   printf("generating single distortion map.  Caution:  This is vastly less efficient than re-using the tpc model once it is set up\n");
0021
0022   bool hasTwin=true; //this flag prompts the code to build both a positive-half and a negative-half for the TPC, reusing as much of the calculations as possible.  It is more efficient to 'twin' one half of the TPC than to recalculate/store the greens functions for both.
0023   TString gainHistName[2]={"hIbfGain_posz","hIbfGFain_negz"};
0024
0025   //and some parameters of the files we're loading:
0026   bool usesChargeDensity=false; //true if source hists contain charge density per bin.  False if hists are charge per bin.
0027   float tpc_chargescale=1.6e-19;//Coulombs per bin unit.
0028   float spacecharge_cm_per_axis_unit=0.1;//cm per histogram axis unit (mm), matching the MDC2 sample from Evgeny.
0029
0030   //file names we'll be filling as we go:
0031   TFile *infile, *gainfile;
0032
0033   TString sourcefilename=inputname;
0034   TString outputfilename=outputname;
0035
0036   //now build the time-consuming part:
0037   AnnularFieldSim *tpc;
0038     tpc=SetupDefaultSphenixTpc(hasTwin,hasSpacecharge, xshift, yshift, zshift);//loads the lookup, fields, etc.
0039
0040   //and the location to plot the fieldslices about:
0041  TVector3 pos=0.5*(tpc->GetOuterEdge()+tpc->GetInnerEdge());;
0042   pos.SetPhi(3.14159);
0043
0044   infile=TFile::Open(sourcefilename.Data(),"READ");
0045
0046
0047   //the total charge is prim + IBF
0048   //if we are doing ADCs, though, we only read the one.
0049   TH3* hCharge=(TH3*)(infile->Get(ibfName));
0050   if (!isAdc){
0051     hCharge->Add((TH3*)(infile->Get(primName)));
0052   }
0053     //hCharge->Scale(70);//Scaleing the histogram spacecharge by 100 times
0054
0055   TString chargestring;
0056
0057   //load the spacecharge into the distortion map generator:
0058   //  void load_spacecharge(TH3F *hist, float zoffset, float chargescale, float cmscale, bool isChargeDensity);
0059   if (!isAdc){
0060     chargestring=Form("%s:(%s+%s)",inputname,ibfName,primName);
0061     tpc->load_spacecharge(hCharge,0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity, chargestring.Data());
0062     if (hasTwin) tpc->twin->load_spacecharge(hCharge,0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity);
0063   }
0064   if (isAdc){ //load digital current using the scaling:
0065     gainfile=TFile::Open(gainName,"READ");
0066     TH2* hGain[2];
0067     hGain[0]=(TH2*)(gainfile->Get(gainHistName[0]));
0068     chargestring=Form("%s:(dc:%s g:%s:%s)",inputname,ibfName,gainName,gainHistName[0].Data());
0069     tpc->load_digital_current(hCharge,hGain[0],tpc_chargescale,spacecharge_cm_per_axis_unit,chargestring.Data());
0070     if (hasTwin) {
0071       hGain[1]=(TH2*)(gainfile->Get(gainHistName[1]));
0072       tpc->twin->load_digital_current(hCharge,hGain[1],tpc_chargescale,spacecharge_cm_per_axis_unit,chargestring.Data());
0073     }
0074   }
0075   //build the electric fieldmap from the chargemap
0076   tpc->populate_fieldmap();
0077   if (hasTwin)  tpc->twin->populate_fieldmap();
0078
0079   //build the distortion maps from the fieldmaps and save it to the output filename.
0080   if (scanSteps){
0081   for(int i=0;i<10;i++){
0082     TString study_filestring=Form("%s.steps%d.hist.root","study_file_changinginterval",50*(i+1));
0083   tpc->GenerateSeparateDistortionMaps(std::string(study_filestring.Data()),50*(i+1),1,1,1,1,true);
0084   //tpc->GenerateSeparateDistortionMaps(outputfilename.Data(),1,1,1,1,false);
0085 }
0086 }
0087   else{
0088    tpc->GenerateSeparateDistortionMaps(std::string(outputfilename.Data()),450,1,1,1,1,false);
0089 }
0090
0091   printf("distortions mapped.\n");
0092   tpc->PlotFieldSlices(outputfilename.Data(),pos, 'E'); //plot the electric field
0093   tpc->PlotFieldSlices(outputfilename.Data(),pos,'B'); //plot the magnetic field
0094   printf("fieldslices plotted.\n");
0095
0096   infile->Close();
0097   printf("input closed.  All done.  Anything else is root's problem.\n");
0098
0099   return;
0100
0101 }
0102
0103
0104 void generate_distortion_map(const char * inputpattern="./evgeny_apr/Smooth*.root", const char *outputfilebase="./apr07_maps/apr07", bool hasSpacecharge=true, bool isDigitalCurrent=false, int nSteps=500){
0105
0106
0107   int maxmaps=10;
0108   int maxmapsperfile=2;
0109
0110
0111   bool hasTwin=true; //this flag prompts the code to build both a positive-half and a negative-half for the TPC, reusing as much of the calculations as possible.  It is more efficient to 'twin' one half of the TPC than to recalculate/store the greens functions for both.
0112   //bool hasSpacecharge=true;
0113
0114   //and some parameters of the files we're loading:
0115   bool usesChargeDensity=false; //true if source hists contain charge density per bin.  False if hists are charge per bin.
0116   float tpc_chargescale=1.6e-19;//Coulombs per bin unit.
0117   float spacecharge_cm_per_axis_unit=0.1;//cm per histogram axis unit, matching the MDC2 sample from Evgeny.
0118
0119   //file names we'll be filling as we go:
0120   TFile *infile;
0121
0122   TString sourcefilename;
0123   TString outputfilename;
0124
0125
0126   //find all files that match the input string (we don't need this yet, but doing it before building the fieldsim saves time if there's an empty directory or something.
0127   TFileCollection *filelist=new TFileCollection();
0128   filelist->Add(inputpattern);
0129   filelist->Print();
0130   printf("Using pattern \"%s\", found %d files to read, eg : %s\n",inputpattern,filelist->GetList()->GetEntries(),((TFileInfo*)(filelist->GetList()->At(0)))->GetCurrentUrl()->GetUrl());//Title());//Print();
0131
0132
0133   // SurveyFiles( filelist);
0134
0135   //now build the time-consuming part:
0136   //AnnularFieldSim *tpc=SetupDefaultSphenixTpc(hasTwin,hasSpacecharge);//loads the lookup, fields, etc.
0137   AnnularFieldSim *tpc;
0138
0139
0140   //previously I flagged digital current and used a different rossegger lookup table.
0141   //now I just resample the histogram in that event.  There is a ~hard problem there to get right:
0142   //in each direction, if the source hist is binned more finely than the dest hist, I need to sum the total charge
0143   //if the source hist is binned more coarsely, I need to interpolate the charge density and multiply by the dest cell volume
0144   //but really, I need to differentiate between the source geometry and field point geometry.  task for another time...
0145   // if (isDigitalCurrent){
0146   //    tpc=SetupDigitalCurrentSphenixTpc(hasTwin,hasSpacecharge);//loads the lookup, fields, etc.
0147   //  }else {
0148     tpc=SetupDefaultSphenixTpc(hasTwin,hasSpacecharge);//loads the lookup, fields, etc.
0149     //  }
0150   //and the location to plot the fieldslices about:
0151  TVector3 pos=0.5*(tpc->GetOuterEdge()+tpc->GetInnerEdge());;
0152   pos.SetPhi(3.14159);
0153
0154
0155
0156
0157
0158   double totalQ=0;
0159
0160   for (int i=0;i<filelist->GetNFiles();i++){
0161    //for each file, find all histograms in that file.
0162     sourcefilename=((TFileInfo*)(filelist->GetList()->At(i)))->GetCurrentUrl()->GetFile();//gross
0163     infile=TFile::Open(sourcefilename.Data(),"READ");
0164     TList *keys=infile->GetListOfKeys();
0165     keys->Print();
0166     int nKeys=infile->GetNkeys();
0167
0168     for (int j=0;j<nKeys;j++){
0169       TObject *tobj=infile->Get(keys->At(j)->GetName());
0170       //if this isn't a 3d histogram, skip it:
0171       bool isHist=tobj->InheritsFrom("TH3");
0172       if (!isHist) continue;
0173       TString objname=tobj->GetName();
0174       if (objname.Contains("IBF")) continue; //this is an IBF-only map we don't want.
0175       //assume this histogram is a charge map.
0176       //load just the averages:
0177       if (hasSpacecharge){
0178     if (isDigitalCurrent){
0179       tpc->load_and_resample_spacecharge(8,36,40,sourcefilename.Data(),tobj->GetName(),0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity);
0180       if (hasTwin) tpc->twin->load_and_resample_spacecharge(8,36,40,sourcefilename.Data(),tobj->GetName(),0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity);
0181
0182     }
0183
0184     //load the spacecharge into the distortion map generator:
0185     tpc->load_spacecharge(sourcefilename.Data(),tobj->GetName(),0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity);
0186     if (hasTwin) tpc->twin->load_spacecharge(sourcefilename.Data(),tobj->GetName(),0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity);
0187
0188     //or load just the average:
0189     //tpc->load_spacecharge(sourcefilename.Data(),"h_Charge_0",0,tpc_chargescale,spacecharge_cm_per_axis_unit, usesChargeDensity);
0190     //printf("Sanity check:  Q has %d elements and dim=%d\n",tpc->q->Length(), tpc->q->dim);
0191     //totalQ=0;
0192     //for (int k=0;k<tpc->q->Length();k++){
0193     //  totalQ+=*(tpc->q->GetFlat(k));
0194     //}
0195     printf("Sanity check:  Total Q in reported region is %E C\n",totalQ);
0196       }
0197       tpc->populate_fieldmap();
0198       if (hasTwin)  tpc->twin->populate_fieldmap();
0199       if (sourcefilename.Contains("verage")){ //this is an IBF map we don't want.
0200     outputfilename=Form("%s.average.%s.%s",outputfilebase,field_string,lookup_string);
0201       } else {
0202     outputfilename=Form("%s.file%d.%s.%s.%s",outputfilebase,i,tobj->GetName(),field_string,lookup_string);
0203       }
0204       printf("%s file has %s hist.  field=%s, lookup=%s. no scaling.\n",
0205          sourcefilename.Data(),tobj->GetName(),field_string,lookup_string);
0206
0207       //TestSpotDistortion(tpc);
0208
0209       //tpc->GenerateDistortionMaps(outputfilename,2,2,2,1,true);
0210
0211
0212       tpc->GenerateSeparateDistortionMaps(std::string(outputfilename.Data()),2,2,2,1,true,nSteps);
0213
0214       printf("distortions mapped.\n");
0215       tpc->PlotFieldSlices(outputfilename.Data(),pos);
0216       tpc->PlotFieldSlices(outputfilename.Data(),pos,'B');
0217       printf("fieldslices plotted.\n");
0218       printf("obj %d: getname: %s  inherits from TH3D:%d \n",j,tobj->GetName(),tobj->InheritsFrom("TH3"));
0219       //break; //rcc temp -- uncomment this to process one hist per file.
0220       if (i>maxmaps) return;
0221     }
0222       infile->Close();
0223   }
0224
0225   return;
0226
0227 }
0228
0229 void TestSpotDistortion(AnnularFieldSim *t){
0230        TVector3 dummy(20.9034,-2.3553,-103.4712);
0231       float dummydest=-103.4752;
0232       t->twin->GetStepDistortion(dummydest,dummy,true,false);
0233       dummy.SetZ(dummy.Z()*-1);
0234       t->GetStepDistortion(-dummydest,dummy,true,false);
0235       return;
0236 }
0237
0238 AnnularFieldSim *SetupDefaultSphenixTpc(bool twinMe, bool useSpacecharge, float xshift, float yshift, float zshift){
0239   //step1:  specify the sPHENIX space charge model parameters
0240   const float tpc_rmin=20.0;
0241   const float tpc_rmax=78.0;
0242   float tpc_deltar=tpc_rmax-tpc_rmin;
0243   const float tpc_z=105.5;
0244   const float tpc_cmVolt=-400*tpc_z; //V =V_CM-V_RDO -- volts per cm times the length of the drift volume.
0245   //const float tpc_magField=0.5;//T -- The old value used in carlos's studies.
0246   //const float tpc_driftVel=4.0*1e6;//cm per s  -- used in carlos's studies
0247   const float tpc_driftVel=8.0*1e6;//cm per s  -- 2019 nominal value
0248   const float tpc_magField=1.4;//T -- 2019 nominal value
0249   const char detgeoname[]="sphenix";
0250
0251
0252   //step 2: specify the parameters of the field simulation.  Larger numbers of
0253   // bins will rapidly increase the memory footprint and compute times.  There
0254   // are some ways to mitigate this by setting a small region of interest, or a
0255   // more parsimonious lookup strategy, specified when AnnularFieldSim() is
0256   // actually constructed below.
0257   int nr=26;//10;//24;//159;//159 nominal
0258   int nr_roi_min=0;
0259   int nr_roi=nr;//10;
0260   int nr_roi_max=nr_roi_min+nr_roi;
0261   int nphi=40;//38;//360;//360 nominal
0262   int nphi_roi_min=0;
0263   int nphi_roi=nphi;//38;
0264   int nphi_roi_max=nphi_roi_min+nphi_roi;
0265   int nz=40;//62;//62 nominal
0266   int nz_roi_min=0;
0267   int nz_roi=nz;
0268   int nz_roi_max=nz_roi_min+nz_roi;
0269
0270   bool realB=true;
0271   bool realE=true;
0272
0273
0274
0275   //step 3:  create the fieldsim object.  different choices of the last few arguments will change how it
0276   // builds the lookup table spatially, and how it loads the spacecharge.  The various start-up steps
0277   // are exposed here so they can be timed in the macro.
0278
0279   AnnularFieldSim *tpc;
0280   if (useSpacecharge){
0281     tpc=    new  AnnularFieldSim(tpc_rmin,tpc_rmax,tpc_z,
0282              nr, nr_roi_min,nr_roi_max,1,2,
0283              nphi,nphi_roi_min, nphi_roi_max,1,2,
0284              nz, nz_roi_min, nz_roi_max,1,2,
0285                  tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::FromFile);
0286   }else{
0287     tpc=    new  AnnularFieldSim(tpc_rmin,tpc_rmax,tpc_z,
0288              nr, nr_roi_min,nr_roi_max,1,2,
0289              nphi,nphi_roi_min, nphi_roi_max,1,2,
0290              nz, nz_roi_min, nz_roi_max,1,2,
0291                  tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::NoSpacecharge);
0292   }
0293
0294   tpc->UpdateEveryN(10);//show reports every 10%.
0295
0296     //load the field maps, either flat or actual maps
0297   tpc->setFlatFields(tpc_magField,tpc_cmVolt/tpc_z);
0298   sprintf(field_string,"flat_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0299
0300   if (realE){
0301     tpc->loadEfield("/sphenix/user/rcorliss/field/externalEfield.ttree.root","fTree");
0302     sprintf(field_string,"realE_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0303   }
0304    if (realB){
0305     tpc->load3dBfield("/sphenix/user/rcorliss/field/sphenix3dmaprhophiz.root","fieldmap",1,-1.4/1.5, xshift, yshift,zshift);
0306         //tpc->loadBfield("sPHENIX.2d.root","fieldmap");
0307     sprintf(field_string,"realB_B%2.1f_E%2.1f_x%2.1f_y%2.1f_z%2.1f",tpc_magField,tpc_cmVolt/tpc_z,xshift, yshift,zshift);
0308   }
0309   if (realE && realB){
0310     sprintf(field_string,"real_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0311   }
0312   printf("set fields.\n");
0313
0314
0315
0316
0317   //load the greens functions:
0318   sprintf(lookup_string,"ross_phi1_%s_phislice_lookup_r%dxp%dxz%d",detgeoname,nr,nphi,nz);
0319   char lookupFilename[300];
0320   //sprintf(lookupFilename,"%s.root",lookup_string);
0321   sprintf(lookupFilename,"/sphenix/user/rcorliss/rossegger/%s.root",lookup_string); //hardcoded for racf
0322   TFile *fileptr=TFile::Open(lookupFilename,"READ");
0323
0324   if (!fileptr){ //generate the lookuptable
0325     //to use the full rossegger terms instead of trivial free-space greens functions, uncomment the line below:
0326     tpc->load_rossegger();
0327     printf("loaded rossegger greens functions.\n");
0328     tpc->populate_lookup();
0329     tpc->save_phislice_lookup(lookupFilename);
0330   } else{ //load it from a file
0331     fileptr->Close();
0332     tpc->load_phislice_lookup(lookupFilename);
0333   }
0334
0335   printf("populated lookup.\n");
0336
0337
0338
0339
0340
0341
0342
0343   //make our twin:
0344   if(twinMe){
0345     AnnularFieldSim *twin;
0346       if (useSpacecharge){
0347     twin=      new  AnnularFieldSim(tpc_rmin,tpc_rmax,-tpc_z,
0348                nr, nr_roi_min,nr_roi_max,1,2,
0349                nphi,nphi_roi_min, nphi_roi_max,1,2,
0350                nz, nz_roi_min, nz_roi_max,1,2,
0351                tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::FromFile);
0352       } else{
0353     twin=      new  AnnularFieldSim(tpc_rmin,tpc_rmax,-tpc_z,
0354                nr, nr_roi_min,nr_roi_max,1,2,
0355                nphi,nphi_roi_min, nphi_roi_max,1,2,
0356                nz, nz_roi_min, nz_roi_max,1,2,
0357                tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::NoSpacecharge);
0358       }    twin->UpdateEveryN(10);//show reports every 10%.
0359
0360     //same magnetic field, opposite electric field
0361     twin->setFlatFields(tpc_magField,-tpc_cmVolt/tpc_z);
0362     //sprintf(field_string,"flat_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0363     //twin->loadBfield("sPHENIX.2d.root","fieldmap");
0364     twin->load3dBfield("/sphenix/user/rcorliss/rossegger/sphenix3dmaprhophiz.root","fieldmap",1,-1.4/1.5, xshift, yshift, zshift);
0365     twin->loadEfield("/sphenix/user/rcorliss/rossegger/externalEfield.ttree.root","fTree",-1);//final '-1' tells it to flip z and the field z coordinate. r coordinate doesn't change, and we assume phi won't either, though the latter is less true.
0366
0367
0368
0369
0370     //borrow the greens functions:
0371     twin->borrow_rossegger(tpc->green,tpc_z);//use the original's green's functions, shift our internal coordinates by tpc_z when querying those functions.
0372     twin->borrow_epartial_from(tpc,tpc_z);//use the original's epartial.  Note that those values ought to be symmetric about z, and since our boundary conditions are translated along with our coordinates, they're completely unchanged.
0373
0374     tpc->set_twin(twin);
0375   }
0376
0377   return tpc;
0378 }
0379
0380
0381 AnnularFieldSim *SetupDigitalCurrentSphenixTpc(bool twinMe, bool useSpacecharge){
0382   //step1:  specify the sPHENIX space charge model parameters
0383   const float tpc_rmin=20.0;
0384   const float tpc_rmax=78.0;
0385   float tpc_deltar=tpc_rmax-tpc_rmin;
0386   const float tpc_z=105.5;
0387   const float tpc_cmVolt=-400*tpc_z; //V =V_CM-V_RDO -- volts per cm times the length of the drift volume.
0388   //const float tpc_magField=0.5;//T -- The old value used in carlos's studies.
0389   //const float tpc_driftVel=4.0*1e6;//cm per s  -- used in carlos's studies
0390   const float tpc_driftVel=8.0*1e6;//cm per s  -- 2019 nominal value
0391   const float tpc_magField=1.4;//T -- 2019 nominal value
0392   const char detgeoname[]="sphenix";
0393
0394   //step 2: specify the parameters of the field simulation.  Larger numbers of
0395   // bins will rapidly increase the memory footprint and compute times.  There
0396   // are some ways to mitigate this by setting a small region of interest, or a
0397   // more parsimonious lookup strategy, specified when AnnularFieldSim() is
0398   // actually constructed below.
0399
0400   int nr=8;//10;//24;//159;//159 nominal
0401   int nr_roi_min=0;
0402   int nr_roi=nr;//10;
0403   int nr_roi_max=nr_roi_min+nr_roi;
0404   int nphi=3*12;//38;//360;//360 nominal
0405   int nphi_roi_min=0;
0406   int nphi_roi=nphi;//38;
0407   int nphi_roi_max=nphi_roi_min+nphi_roi;
0408   int nz=40;//62;//62 nominal
0409   int nz_roi_min=0;
0410   int nz_roi=nz;
0411   int nz_roi_max=nz_roi_min+nz_roi;
0412
0413   bool realB=true;
0414   bool realE=true;
0415
0416
0417   //step 3:  create the fieldsim object.  different choices of the last few arguments will change how it
0418   // builds the lookup table spatially, and how it loads the spacecharge.  The various start-up steps
0419   // are exposed here so they can be timed in the macro.
0420   AnnularFieldSim *tpc;
0421   if (useSpacecharge){
0422     tpc=    new  AnnularFieldSim(tpc_rmin,tpc_rmax,tpc_z,
0423              nr, nr_roi_min,nr_roi_max,1,2,
0424              nphi,nphi_roi_min, nphi_roi_max,1,2,
0425              nz, nz_roi_min, nz_roi_max,1,2,
0426                  tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::FromFile);
0427   }else{
0428     tpc=    new  AnnularFieldSim(tpc_rmin,tpc_rmax,tpc_z,
0429              nr, nr_roi_min,nr_roi_max,1,2,
0430              nphi,nphi_roi_min, nphi_roi_max,1,2,
0431              nz, nz_roi_min, nz_roi_max,1,2,
0432                  tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::NoSpacecharge);
0433   }
0434
0435   tpc->UpdateEveryN(10);//show reports every 10%.
0436
0437     //load the field maps, either flat or actual maps
0438   tpc->setFlatFields(tpc_magField,tpc_cmVolt/tpc_z);
0439   sprintf(field_string,"flat_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0440
0441   if (realE){
0442     tpc->loadEfield("externalEfield.ttree.root","fTree");
0443     sprintf(field_string,"realE_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0444   }
0445    if (realB){
0446     tpc->load3dBfield("sphenix3dmaprhophiz.root","fieldmap",1,-1.4/1.5);
0447         //tpc->loadBfield("sPHENIX.2d.root","fieldmap");
0448     sprintf(field_string,"realB_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0449   }
0450   if (realE && realB){
0451     sprintf(field_string,"real_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0452   }
0453   printf("set fields.\n");
0454
0455
0456
0457
0458   //load the greens functions:
0459   sprintf(lookup_string,"ross_phi1_%s_phislice_lookup_r%dxp%dxz%d",detgeoname,nr,nphi,nz);
0460   char lookupFilename[200];
0461   sprintf(lookupFilename,"%s.root",lookup_string);
0462   TFile *fileptr=TFile::Open(lookupFilename,"READ");
0463
0464   if (!fileptr){ //generate the lookuptable
0465     //to use the full rossegger terms instead of trivial free-space greens functions, uncomment the line below:
0466     tpc->load_rossegger();
0467     printf("loaded rossegger greens functions.\n");
0468     tpc->populate_lookup();
0469     tpc->save_phislice_lookup(lookupFilename);
0470   } else{ //load it from a file
0471     fileptr->Close();
0472     tpc->load_phislice_lookup(lookupFilename);
0473   }
0474
0475   printf("populated lookup.\n");
0476
0477
0478
0479
0480
0481
0482
0483   //make our twin:
0484   if(twinMe){
0485     AnnularFieldSim *twin;
0486       if (useSpacecharge){
0487     twin=      new  AnnularFieldSim(tpc_rmin,tpc_rmax,-tpc_z,
0488                nr, nr_roi_min,nr_roi_max,1,2,
0489                nphi,nphi_roi_min, nphi_roi_max,1,2,
0490                nz, nz_roi_min, nz_roi_max,1,2,
0491                tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::FromFile);
0492       } else{
0493     twin=      new  AnnularFieldSim(tpc_rmin,tpc_rmax,-tpc_z,
0494                nr, nr_roi_min,nr_roi_max,1,2,
0495                nphi,nphi_roi_min, nphi_roi_max,1,2,
0496                nz, nz_roi_min, nz_roi_max,1,2,
0497                tpc_driftVel, AnnularFieldSim::PhiSlice, AnnularFieldSim::NoSpacecharge);
0498       }    twin->UpdateEveryN(10);//show reports every 10%.
0499
0500     //same magnetic field, opposite electric field
0501     twin->setFlatFields(tpc_magField,-tpc_cmVolt/tpc_z);
0502     //sprintf(field_string,"flat_B%2.1f_E%2.1f",tpc_magField,tpc_cmVolt/tpc_z);
0503     //twin->loadBfield("sPHENIX.2d.root","fieldmap");
0504     twin->load3dBfield("sphenix3dmaprhophiz.root","fieldmap",1,-1.4/1.5);
0505     twin->loadEfield("externalEfield.ttree.root","fTree",-1);//final '-1' tells it to flip z and the field z coordinate. r coordinate doesn't change, and we assume phi won't either, though the latter is less true.
0506
0507
0508
0509
0510     //borrow the greens functions:
0511     twin->borrow_rossegger(tpc->green,tpc_z);//use the original's green's functions, shift our internal coordinates by tpc_z when querying those functions.
0512     twin->borrow_epartial_from(tpc,tpc_z);//use the original's epartial.  Note that those values ought to be symmetric about z, and since our boundary conditions are translated along with our coordinates, they're completely unchanged.
0513
0514     tpc->set_twin(twin);
0515   }
0516
0517   return tpc;
0518 }
0519
0520 void  SurveyFiles(TFileCollection *filelist){
0521    TFile *infile;
0522
0523   TString sourcefilename;
0524   TString outputfilename;
0525  //run a check of the files we'll be looking at:
0526   float integral_sum=0;
0527   int nhists=0;
0528   for (int i=0;i<filelist->GetNFiles();i++){
0529    //for each file, find all histograms in that file.
0530     sourcefilename=((TFileInfo*)(filelist->GetList()->At(i)))->GetCurrentUrl()->GetFile();//gross
0531     printf("file %d: %s\n", i, sourcefilename.Data());
0532     infile=TFile::Open(sourcefilename.Data(),"READ");
0533     TList *keys=infile->GetListOfKeys();
0534     //keys->Print();
0535     int nKeys=infile->GetNkeys();
0536
0537     for (int j=0;j<nKeys;j++){
0538       TObject *tobj=infile->Get(keys->At(j)->GetName());
0539       //if this isn't a 3d histogram, skip it:
0540       bool isHist=tobj->InheritsFrom("TH3");
0541       if (!isHist) continue;
0542       TString objname=tobj->GetName();
0543       printf(" hist %s ",objname.Data());
0544       if (objname.Contains("IBF")) {
0545     printf(" is IBF only.\n");
0546     continue; //this is an IBF map we don't want.
0547       }
0548       float integral=((TH3D*)tobj)->Integral();
0549       integral_sum+=integral;
0550       nhists+=1;
0551     printf(" will be used.  Total Q=%3.3E (ave=%3.3E)\n",integral,integral_sum/nhists);
0552       //assume this histogram is a charge map.
0553       //load just the averages:
0554
0555       //break; //rcc temp -- uncomment this to process one hist per file.
0556       //if (i>maxmaps) return;
0557     }
0558       infile->Close();
0559   }
0560   return;
0561 }