scansObsids = [1342202254]
cscansObsids = [1342190616, 1342202090]

galactic = True
debug = False
deglitch = True
nSigmaDeglitch = 3

"""
########################### SECTION 1 ##########################################
################################################################################
########################### PROCESSING #########################################
"""
print " -- Executing JScanam -- " 

print " Processing the following obsid pair: " + str(scansObsids[0]) + " " + str(cscansObsids[0])
   
"""
 Load Level 1 data from the HSA
"""
print " Loading Level 1 data "

for i in range(len(scansObsids)):
   scansObs = getObservation(scansObsids[i], useHsa=True, instrument="PACS")
   level1 = PacsContext(scansObs.level1)
   if i == 0:
      scans = level1.averaged.getCamera(camera).product.selectAll() 
   else:
      scans.join(level1.averaged.getCamera(camera).product.selectAll()) 
   blueFilter1 = scansObs.meta["blue"].value

for i in range(len(cscansObsids)):
   cscansObs = getObservation(cscansObsids[i], useHsa=True, instrument="PACS")
   level1 = PacsContext(cscansObs.level1)
   if i == 0:
      cscans = level1.averaged.getCamera(camera).product.selectAll() 
   else:
      cscans.join(level1.averaged.getCamera(camera).product.selectAll()) 
   blueFilter2 = cscansObs.meta["blue"].value

calTree = getCalTree(obs=cscansObs)

if (camera == "blue" and blueFilter1 != blueFilter2):
   print ""
   print "  ------------------------------------------------------------------"
   print "  ALERT!!! You are trying to combine two blue observations obtained "  
   print "  with different filters, which is most probably wrong!             "
   print "  ------------------------------------------------------------------"
   print ""

"""
 Remove turn arounds
"""
print " Removing turn arounds "
scans = scanamorphosRemoveTurnarounds(scans, limit=50.0, debug=debug)
cscans = scanamorphosRemoveTurnarounds(cscans, limit=50.0, debug=debug)

"""
 Mask long term glitches. 
 This task produces a mask called Scanam_LongTermGlitchMask. You should check 
 this mask and if the results are not optimal (some glitches are not detected 
 or some sources are flagged as glitches), you can try to modify the parameter 
 stepAfter in order to get better results.
"""
print " Masking long term glitches "
scans = scanamorphosMaskLongTermGlitches(scans, stepAfter=20, galactic=galactic, calTree=calTree, debug=debug)
cscans = scanamorphosMaskLongTermGlitches(cscans, stepAfter=20, galactic=galactic, calTree=calTree, debug=debug)

"""
 Save the scans and cross scans for later use
"""
print " Saving a copy of the scans and cross scans "
scansCopy = scans.copy()
cscansCopy = cscans.copy()

"""
 Subtract the baseline of every scanleg in every pixel
"""
print " Initial baseline subtraction / per scanleg "
scans = scanamorphosScanlegBaselineFitPerPixel(scans, nSigma=2)
cscans = scanamorphosScanlegBaselineFitPerPixel(cscans, nSigma=2)

"""
 Create the source mask.
 Modify nSigma until you get an optimal mask. The mask should cover only a small
 fraction of the map (<~30%). It's not necessary that all the faint emission is
 masked, only the brightest regions.
"""
print " Creating a source mask "
scans.join(cscans)
del(cscans)
sourceImage, scans = scanamorphosCreateSourceMask(scans, nSigma=4.0, createMask=False, galactic=galactic, calTree=calTree, debug=debug)

"""
 Replace the scans and cross scans by the saved ones
"""
print " Loading the saved scans and cross scans from the temporal pool "
del(scans)
scans = scansCopy
cscans = cscansCopy
del(scansCopy, cscansCopy)
System.gc()
 
"""
 Add the mask to the saved scans and cross scans
"""
print " Adding the source mask to the scans and cross scans "
maskImage, scans = scanamorphosCreateSourceMask(scans, inputImage=sourceImage, createMask=True, calTree=calTree, debug=debug)
maskImage, cscans = scanamorphosCreateSourceMask(cscans, inputImage=sourceImage, createMask=True, calTree=calTree, debug=debug)

"""
 Baseline subtraction. 
 Here we use galactic=True, because we only want to remove an offset
"""
print " Baseline subtraction "
scans = scanamorphosBaselineSubtraction(scans, galactic=True, debug=debug)
cscans = scanamorphosBaselineSubtraction(cscans, galactic=True, debug=debug)

"""
 Baseline pre-processing
"""
print " Baseline pre-processing "
scans = scanamorphosBaselinePreprocessing(scans, debug=debug)
cscans = scanamorphosBaselinePreprocessing(cscans, debug=debug)

"""
 Baseline subtraction
"""
print " Baseline subtraction "
scans = scanamorphosBaselineSubtraction(scans, galactic=galactic, debug=debug)
cscans = scanamorphosBaselineSubtraction(cscans, galactic=galactic, debug=debug)

"""
 Destriping
"""
print " Destriping "
scans, cscans = scanamorphosDestriping(scans, cscans, iterations=6, calTree=calTree, debug=debug)

"""
 Align the frames to have the same background and create the no-ovelap mask
"""
print " Aligning the scans and cross scans background reference "
scans, cscans = scanamorphosAlignFrames(scans, cscans, calTree=calTree)

"""
 Merge the scans and cross scans
"""
print " Merging scans and cross scans in a frame product "
System.gc()
mergedScans = scans
mergedScans.join(cscans)
del(scans, cscans)
System.gc()

"""
 Deglitch the merged scans
 
 We use nSigma=5, and not a lower value, to avoid masking signal with strong drifts  
 that will be corrected by the scanamorphosIndividualDrifts task
"""
print " Deglitching the merged scans "
scanamorphosDeglitch(mergedScans, nSigma=5, calTree=calTree, debug=debug)

"""
 Remove individual drifts
"""
print " Removing individual drifts "
mergedScans = scanamorphosIndividualDrifts(mergedScans, calTree=calTree, debug=debug)

"""
 Deglitch the merged scans again with the user provided nSigmaDeglitch

 Be careful setting nSigmaDeglitch to very small values, because it can
 mask bright sources. You should always check the mask called Scanamorphos_GlitchMask
 to make sure that the task is masking only real glitches.
"""
if(deglitch):
   print " Deglitching the merged scans "
   scanamorphosDeglitch(mergedScans, nSigma=nSigmaDeglitch, calTree=calTree, debug=debug)

"""
  Write the fits files
"""
print " Saving the fits file "
cameraName = camera
if(camera == "blue"):
   if(blueFilter1 == "blue2"):
      cameraName = "green"

outputFitsFile = dataDir + str(scansObsids[0]) +  "-" + str(cscansObsids[0])  +  "-" + cameraName + ".fits"
FitsArchive().save(outputFitsFile, mergedScans)

del(mergedScans)
System.gc()